Changeset 1988 for trunk

Timestamp:

Aug 28, 2018, 4:22:24 PM (6 years ago)

Author:

jleconte

Message:

28/08/2018 == JL

We now shift the radiative model top from p=0 to the middle of the last physical layer. This is done by changing pmid and plevrad in callcorrk and some corrections need to be done in gfluxv.
This seems to get rid of the aratic temperature behavior in the last two layers of the model (especially on the night side on synchronous planets).
Additional speedup corrections have been made in gfluxi that change nothing to the result.
Finally, if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv and changes in last commit.)
This has been done for water ice in aeropacity, but same correction should probably be done for other aerosol types.

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/phystd/aeropacity.F90 (modified) (3 diffs)
• libf/phystd/callcorrk.F90 (modified) (3 diffs)
• libf/phystd/gfluxi.F (modified) (2 diffs)
• libf/phystd/gfluxv.F (modified) (3 diffs)

trunk/LMDZ.GENERIC/README

@@ -1378 +1378 @@
 Tauaero and tauray are set to 0 (a small value for rayleigh because the code crashes otherwise) in the 4 first semilayers in optcv, but not optci. 
 This solves random variations of the sw heating at the model top. 
+
+== 28/08/2018 == JL
+We now shift the radiative model top from p=0 to the middle of the last physical layer. This is done by changing pmid and plevrad in callcorrk and some corrections need to be done in gfluxv. 
+This seems to get rid of the aratic temperature behavior in the last two layers of the model (especially on the night side on synchronous planets).
+Additional speedup corrections have been made in gfluxi that change nothing to the result.
+Finally, if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv and changes in last commit.)
+This has been done for water ice in aeropacity, but same correction should probably be done for other aerosol types.
+

trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90

@@ -219 +219 @@
 
                do ig=1, ngrid
-                  do l=1,nlayer-1 ! to stop the rad tran bug
-
+                  !do l=1,nlayer-1 ! to stop the rad tran bug
+                  do l=1,nlayer !JL18 if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv)
+                                ! same correction should b-probably be done for other aerosol types.
                      aerosol(ig,l,iaer) =                                    & !modification by BC
                           (  0.75 * QREFvis3d(ig,l,iaer) /        &
@@ -236 +237 @@
                   call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
                   do ig=1, ngrid
-                     do l=1,nlayer-1 ! to stop the rad tran bug
+                     !do l=1,nlayer-1 ! to stop the rad tran bug
+                     do l=1,nlayer !JL18 if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv)
                         CLFtot  = max(totcloudfrac(ig),0.01)
                         aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
@@ -244 +246 @@
                else
                   do ig=1, ngrid
-                     do l=1,nlayer-1 ! to stop the rad tran bug
+                     !do l=1,nlayer-1 ! to stop the rad tran bug
+                     do l=1,nlayer !JL18 if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv)
                         CLFtot  = CLFfixval
                         aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot

trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

@@ -486 +486 @@
             tauaero(1,iaer)          = tauaero(2,iaer)
             !tauaero(1,iaer)          = 0.
-            
+            !JL18 at time of testing, the two above conditions gave the same results bit for bit. 
+            
          end do ! naerkind
 
@@ -665 +666 @@
       
       plevrad(1) = 0.
-      plevrad(2) = 0.   !! Trick to have correct calculations of fluxes in gflux(i/v).F, but the pmid levels are not impacted by this change.
+!      plevrad(2) = 0.   !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all. 
 
       tlevrad(1) = tlevrad(2)
       tlevrad(2*nlayer+1)=tsurf(ig)
       
-      pmid(1) = max(pgasmin,0.0001*plevrad(3))
+      pmid(1) = pplay(ig,nlayer)/scalep
       pmid(2) =  pmid(1)
 
@@ -899 +900 @@
 !     These are values at top of atmosphere
          dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
-             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1)))
+             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
          dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
-             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1)))
+             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
 
 

trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F

@@ -117 +117 @@
         TERM    = UBARI/(1.0D0-W0(L)*COSBAR(L))
 
-C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
-C       BOTTOM OF THE LAYER.  THAT IS AT DTAU OPTICAL DEPTH
-
-        CP(L) = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM
-        CM(L) = B0(L)+B1(L)*DTAU(L) -B1(L)*TERM
 
 C       CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
@@ -128 +123 @@
         CPM1(L) = B0(L)+B1(L)*TERM
         CMM1(L) = B0(L)-B1(L)*TERM
+
+C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
+C       BOTTOM OF THE LAYER.  THAT IS AT DTAU OPTICAL DEPTH.
+C       JL18 put CP and CM after the calculation of CPM1 and CMM1 to avoid unecessary calculations. 
+
+        CP(L) = CPM1(L) +B1(L)*DTAU(L) 
+        CM(L) = CMM1(L) +B1(L)*DTAU(L) 
       END DO
  

trunk/LMDZ.GENERIC/libf/phystd/gfluxv.F

@@ -240 +240 @@
 C     FLUX AT THE Ptop layer
 
-      EP    = 1.0
-      EM    = 1.0
+!      EP    = 1.0
+!      EM    = 1.0
+C JL18 correction to account for the fact that the radiative top is not at zero optical depth.
+      EXPTRM(L) = MIN(TAUMAX,LAMDA(L)*(TAUCUM(2)))
+      EP = EXP(EXPTRM(L))
+      EM    = 1.0/EP
       G4    = 1.0-G3(1)
       DENOM = LAMDA(1)**2 - 1./UBAR0**2
@@ -260 +264 @@
 C     AT THE MIDDLE OF THE LAYER.
 
-      CPMID  = AP
-      CMMID  = AM
+C      CPMID  = AP
+C      CMMID  = AM
+C JL18 correction to account for the fact that the radiative top is not at zero optical depth.
+      TAUMID   = TAUCUM(2)
+      CPMID = AP*EXP(-TAUMID/UBAR0)
+      CMMID = AM*EXP(-TAUMID/UBAR0)
 
       FLUXUP = XK1(1)*EP + GAMA(1)*XK2(1)*EM + CPMID
@@ -268 +276 @@
 C     ADD THE DIRECT FLUX TO THE DOWNWELLING TERM
 
-      fluxdn = fluxdn+UBAR0*F0PI*EXP(-MIN(TAUCUM(1),TAUMAX)/UBAR0)
+!      fluxdn = fluxdn+UBAR0*F0PI*EXP(-MIN(TAUCUM(1),TAUMAX)/UBAR0)
+!JL18 the line above assumed that the top of the radiative model was P=0
+!   it seems to be better for the IR to use the middle of the last physical layer as the radiative top. 
+!   so we correct the downwelling flux below for the calculation of the heating rate
+      fluxdn = fluxdn+UBAR0*F0PI*EXP(-TAUCUM(2)/UBAR0)
 
 C     This is for the "special" bottom layer, where we take