From 148f9bccee18397142cb12ca68f9bf62cdf39cee Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Wed, 11 May 2022 16:51:23 +0100 Subject: [PATCH] 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 | 2 +- exoplasim/plasim/src/aerocore.f90 | 12 ++++++------ exoplasim/plasim/src/plasim.f90 | 18 ++++++++---------- exoplasim/pyburn.py | 3 ++- 4 files changed, 17 insertions(+), 18 deletions(-) diff --git a/exoplasim/plasim/run/example.nl b/exoplasim/plasim/run/example.nl index c0f7462..4e8cef4 100644 --- a/exoplasim/plasim/run/example.nl +++ b/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 diff --git a/exoplasim/plasim/src/aerocore.f90 b/exoplasim/plasim/src/aerocore.f90 index c10aa7d..69679ce 100755 --- a/exoplasim/plasim/src/aerocore.f90 +++ b/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 diff --git a/exoplasim/plasim/src/plasim.f90 b/exoplasim/plasim/src/plasim.f90 index a0796bb..103959a 100644 --- a/exoplasim/plasim/src/plasim.f90 +++ b/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 diff --git a/exoplasim/pyburn.py b/exoplasim/pyburn.py index 5daed67..b1c49ea 100644 --- a/exoplasim/pyburn.py +++ b/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 = {}