Now running on one PRO without numerical instabilty and output container showing up in .nc file. Some negatives in mmr and cube in .nc is empty. Still debugging.

exoplasim/plasim/run/example.nl

@@ -1,4 +1,4 @@
-code=129,130,131,132,133,134,135,137,138,139,140,141,142,143,144,146,147,148,149,151,152,155,157,160,161,162,164,167,169,172,175,176,177,178,179,182,183,173,203,204,205,207,210,211,218,221,230,232,259,260,261,262,263,264,266,267,301,302,318,319,320,321,322,323,324,325,326,327,328,329
+code=129,130,131,132,133,134,135,137,138,139,140,141,142,143,144,146,147,148,149,151,152,155,157,160,161,162,164,167,169,172,175,176,177,178,179,182,183,173,203,204,205,207,210,211,218,221,230,232,259,260,261,262,263,264,266,267,301,302,318,319,320,321,322,323,324,325,326,327,328,329,410
 vtype=s
 htype=g
 mean=3

exoplasim/plasim/src/aerocore.f90

@@ -693,7 +693,7 @@
 ! Compute fluxes in the vertical direction
       call FZPPM(qz,fz,IM,JM,NL,daero,W,delp,KORD)
 !****6***0*********0*********0*********0*********0*********0**********72
-
+      write(nud,*) 'Max and min mmr =',maxval(mmr), minval(mmr)
 !****6***0*********0*********0*********0*********0*********0**********72
 ! Compute vertical flux due to gravitational settling 
       call gsettle(qz,im,jm,nl,rhog,vels,nud,gz)
@@ -737,13 +737,13 @@
 !MIC$ do all autoscope
 !MIC$* private(i,j,k,sum1,sum2)
 
-! For k=1, only subtract the flux due to downward settling (plus all the advection fluxes)
+! For k=1, update top level by subtracting downward flux from this level 
 	  do j=j1,j2 
 		do i=1,IM
 			daero(i,j,1) = daero(i,j,1) +  fx(i,j,1) - fx(i+1,j,1)                   &
                             + (fy(i,j,1) - fy(i,j+1,1))*acosp(j)/dlat(j) &
-                            +  fz(i,j,1) - fz(i,j,2) &
-							- gz(i,j,1)*ga
+                            +  fz(i,j,1) - fz(i,j,2)	 &
+                            - gz(i,j,1)*ga							
 		enddo ! Lon loop
 	  enddo ! Lat loop
 	  
@@ -916,8 +916,8 @@
 	rt_mair = SQRT(mair)
 	rt_pi = SQRT(PI)
 	rt_kb = SQRT(kb)
-	sq_dair = dair**2
-	temp_eps = (temp/eps)*(4/25)
+	sq_dair = dair*dair
+	temp_eps = (temp/eps)**(4/25)
 	coeff = (5./16.)*(1./1.22)*(1./PI)*rt_pi*rt_mair*rt_kb/sq_dair
 	
 	mu = coeff*rt_temp*temp_eps

exoplasim/plasim/src/plasim.f90

@@ -3088,7 +3088,7 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
       dvdt(:,:)=0.
       dtdt(:,:)=0.
       dqdt(:,:)=0.
-!	  mmrt(:,:)=0.
+	  mmrt(:,:)=0.
 
 !
 !     transform to gridpoint domain
@@ -3120,6 +3120,8 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
 !
       if (nsela > 0 .and. nkits == 0) then	  
        if (nqspec == 0) then
+	    write(nud,*) 'Semi-Lagrangian q running'
+	    flush(nud)
         dqt(:,:) = dq(:,:) ! Save old value of q
         call mpgagp(zgq,dq,NLEV)
         if (mypid == NROOT) then
@@ -3278,26 +3280,22 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
 	  ! endif
 
 	  
-      if (nsela > 0 .and. nkits == 0) then
-	   if (NAERO > 0) then
-        mmrt(:,:) = mmr(:,:) ! Save old value of q
+      if (nsela == 1 .and. nkits == 0 .and. NAERO > 0) then
+        mmrt(:,:) = mmr(:,:) ! Save old value of mmr
         call mpgagp(zmmr,mmr,NLEV)
         if (mypid == NROOT) then
          do jlat = 1 , NLAT
           daeros(:,NLAT+1-jlat,:,1) = zmmr(:,jlat,:)
          enddo ! jlat
         endif ! mypid
-       endif ! naero
-       call aero_main
-       if (NAERO > 0) then
+        call aero_main
         if (mypid == NROOT) then
          do jlat = 1 , NLAT
           zmmr(:,jlat,:) = daeros(:,NLAT+1-jlat,:,1)
          enddo
         endif ! mypid
         call mpscgp(zmmr,mmr,NLEV)
-        mmrt(:,:) = (mmr(:,:) - mmrt(:,:)) / deltsec !  q advection term
-       endif ! nqspec	  
+        mmrt(:,:) = (mmr(:,:) - mmrt(:,:)) / deltsec !  q advection term       endif ! nqspec	  
       endif ! nkits 
 
 !
@@ -3393,7 +3391,7 @@ plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
       call mpsumsc(szf,szt,NLEV)
       if (nqspec == 1) call mpsumsc(sqf,sqt,NLEV)
       if (nqspec == 0) dq(:,:) = dq(:,:) + dqdt(:,:) * deltsec
-!	  if (nsela == 1 .and. NAERO > 0) mmr(:,:) = mmr(:,:) + mmrt(:,:) * deltsec
+	  if (nsela == 1 .and. NAERO > 0) mmr(:,:) = mmr(:,:) + mmrt(:,:) * deltsec
 
       return
       end

exoplasim/pyburn.py

@@ -153,7 +153,8 @@ ilibrary = {"110":["mld"  ,"mixed_layer_depth"               ,"m"          ],
             "407":["dftd"  ,"longwave_down"                  ,"W m-2"      ],
             "408":["dtdt"  ,"rad_heating_rate"               ,"K s-1"      ],
             "409":["dfdz"  ,"flux_convergence"               ,"W m-3"      ],
-			"410":["mmr"   ,"aerosol_mass_mixing_ratio"      ,"kg kg-1"    ]}
+			"410":["mmr"   ,"aerosol_mass_mixing_ratio"      ,"kg kg-1"    ],
+			"411":["mmrt"  ,"aerosol_mmr_tendency"           ,"kg kg-1 s-1"]}
 
 #dictionary that can be searched by string codes
 slibrary = {}