Changeset 4146 for trunk/LMDZ.GENERIC/libf/phygeneric/turbdiff_mod.F90

Timestamp:

Mar 19, 2026, 2:35:46 PM (10 days ago)

Author:

gmilcareck

Message:

Thermodynamics update on LMDZ.GENERIC

File:

• trunk/LMDZ.GENERIC/libf/phygeneric/turbdiff_mod.F90 (modified) (12 diffs)

trunk/LMDZ.GENERIC/libf/phygeneric/turbdiff_mod.F90

@@ -8 +8 @@
           ptimestep,pcapcal,                    &   
           pplay,pplev,pzlay,pzlev,pz0,                 &
-          pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf,       &
+          pu,pv,pt,ph,ppopsk,pq,ptsrf,pemis,pqsurf,       &
           pdtfi,pdqfi,pfluxsrf,            &
           Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2,  &
@@ -17 +17 @@
       use surfdat_h, only: dryness
       use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
-      use comcstfi_mod, only: rcp, g, r, cpp
+      use comcstfi_mod, only: g, rd_ref, cppd_ref
       use callkeys_mod, only: water,tracer,nosurf,kmixmin
       use turb_mod, only : ustar
@@ -61 +61 @@
       REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
       REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
-      REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
+      REAL,INTENT(IN) :: pt(ngrid,nlay),ph(ngrid,nlay),ppopsk(ngrid,nlay)
       REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
       REAL,INTENT(IN) :: pemis(ngrid)
@@ -173 +173 @@
          DO ig=1,ngrid
             zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
-            zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
+            zExner(ig,ilay)=ppopsk(ig,ilay)
             zovExner(ig,ilay)=1./ppopsk(ig,ilay)
          ENDDO
       ENDDO
 
-      zcst1=4.*g*ptimestep/(R*R)
+      zcst1=4.*g*ptimestep/(rd_ref*rd_ref)
       DO ilev=2,nlev-1
          DO ig=1,ngrid
@@ -186 +186 @@
       ENDDO
       DO ig=1,ngrid
-         zb0(ig,1)=ptimestep*pplev(ig,1)/(R*ptsrf(ig))
+         zb0(ig,1)=ptimestep*pplev(ig,1)/(rd_ref*ptsrf(ig))
       ENDDO
       dqsdif_total(:)=0.0
@@ -199 +199 @@
             zv(ig,ilev)=pv(ig,ilev)
             zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
-            zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
+            zh(ig,ilev)=ph(ig,ilev)!+pdtfi(ig,ilev)*ptimestep*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
         ENDDO
       ENDDO
@@ -248 +248 @@
 !     ------------------------------------------------------ 
 
-      call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay,pu,pv,pq,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature 
-      
+      call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay,pplev,pplay, & 
+                   pu,pv,zt,pq,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated
+                       
 !     Adding eddy mixing to mimic 3D general circulation in 1D
 !     R. Wordsworth & F. Forget (2010)
@@ -273 +274 @@
 !               zkv(ig,ilev) = max(kzz_eddy(ilev),zkv(ig,ilev)) 
 !            end do
-
             do ig=1,ngrid
                zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
@@ -436 +436 @@
 
          DO ig=1,ngrid
-            z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) &
+            z1(ig)=pcapcal(ig)*ptsrf(ig)+cppd_ref*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) &
                 + pfluxsrf(ig)*ptimestep + zdplanck(ig)*ptsrf(ig) 
-            z2(ig) = pcapcal(ig)+zdplanck(ig)+cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) 
+            z2(ig) = pcapcal(ig)+zdplanck(ig)+cppd_ref*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) 
             ztsrf(ig) = z1(ig) / z2(ig)
             pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
@@ -573 +573 @@
                   zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
 
-                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1)   &
+                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cppd_ref*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1)   &
                       + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep                       &
                       + zfluxq(ig,1)*dryness(ig)*RLVTT*((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
 
-                  z2(ig) = pcapcal(ig) + cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
+                  z2(ig) = pcapcal(ig) + cppd_ref*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
                       + zdplanck(ig)+zfluxq(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*(1.0-zdq(ig,1))
 
@@ -598 +598 @@
 
                       !recompute surface temperature  
-                      z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1)   &
+                      z1(ig) = pcapcal(ig)*ptsrf(ig) +cppd_ref*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1)   &
                         + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep                       &
                         + RLVTT*dqsdif_total(ig)
-                      z2(ig) = pcapcal(ig) + cpp*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
+                      z2(ig) = pcapcal(ig) + cppd_ref*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
                         + zdplanck(ig)
                       ztsrf(ig) = z1(ig) / z2(ig)
@@ -720 +720 @@
       
       DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
-         sensibFlux(ig)=cpp*zfluxt(ig,1)/ptimestep*(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
+         sensibFlux(ig)=cppd_ref*zfluxt(ig,1)/ptimestep*(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
       ENDDO