Changeset 2097 for trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

Timestamp:

Feb 11, 2019, 5:53:10 PM (6 years ago)

Author:

jvatant

Message:

+ Fix some ambiguities about effective gravity field. Choice is made to keep coherence between dynamics and physics by default,

but having correct altitudes for chemistry - to have correct fluxes - no matter what's in physics.
Then chemistry is "recast" on the dynphy grid.

+ Also correct the fact that chemistry requires altitutde with ref to the geoid not to the local surface -> pphi+phhis in the calculation.
+ Both points above lead to new zlay_eff/zlev_eff altitudes variables in physiq_mod sent to chemistry.

We still need to think about it for microphysics.

--JVO

File:

• trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (modified) (10 diffs)

trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

@@ -260 +260 @@
       real zlss                      ! Sub solar point longitude (radians).
       real zday                      ! Date (time since Ls=0, calculated in sols).
-      real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers.
-      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries.
-
+      real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers (ref : local surf).
+      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries (ref : local surf).
+      real zzlay_eff(ngrid,nlayer)   ! Effective altitude at the middle of the atmospheric layers (ref : geoid ).
+      real zzlev_eff(ngrid,nlayer+1) ! Effective altitude at the atmospheric layer boundaries ( ref : geoid ).
+      
 ! TENDENCIES due to various processes :
 
@@ -704 +706 @@
       
         do l=1,nlayer         
-           zzlay(:,l) = pphi(:,l) / gzlat(:,l)
+           zzlay(:,l) = pphi(:,l) / gzlat(:,l) ! Reference = local surface
         enddo
         
@@ -710 +712 @@
       
         do l=1,nlayer
-           zzlay(:,l) = g*rad*rad / ( g*rad - pphi(:,l) ) - rad
+           zzlay(:,l) = g*rad*rad / ( g*rad - ( pphi(:,l) + phisfi(:) )) - rad
            gzlat(:,l) = g*rad*rad / ( rad + zzlay(:,l) )**2
         end do
@@ -718 +720 @@
       zzlev(:,1)=0.
       zzlev(:,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
+      ! JVO 19 : This altitude is indeed dummy for the GCM and fits ptop=0
+      ! but for upper chemistry that's a pb -> we anyway redefine it just after ..
 
       do l=2,nlayer
@@ -727 +731 @@
       enddo     
 
-      ! Zonal averages needed for chemistry
-      ! ~~~~~~~ Taken from old Titan ~~~~~~
-      if (moyzon_ch) then
-         
-         if (eff_gz .eqv. .false.) then
-           zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
-         else ! if we take into account effective altitude-dependent gravity field
-           zzlaybar(1,:)=g*rad*rad/(g*rad-(zphibar(1,:)+zphisbar(1)))-rad
-         endif
+      ! Effective altitudes ( eg needed for chemistry ) with correct g, and with reference to the geoid
+      ! JVO 19 : We shall always have correct altitudes in chemistry no matter what's in physics
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (moyzon_ch) then ! Zonal averages
+         
+         zzlaybar(1,:)=g*rad*rad/(g*rad-(zphibar(1,:)+zphisbar(1)))-rad   ! reference = geoid 
          zzlevbar(1,1)=zphisbar(1)/g       
          DO l=2,nlayer
@@ -747 +748 @@
             if (latitude(ig).ne.latitude(ig-1)) then        
                DO l=1,nlayer   
-                  if (eff_gz .eqv. .false.) then
-                    zzlaybar(ig,l)=(zphibar(ig,l)+zphisbar(ig))/g
-                  else ! if we take into account effective altitude-dependent gravity field
-                    zzlaybar(ig,l)=g*rad*rad/(g*rad-(zphibar(ig,l)+zphisbar(ig)))-rad
-                  endif            
+                  zzlaybar(ig,l)=g*rad*rad/(g*rad-(zphibar(ig,l)+zphisbar(ig)))-rad
                ENDDO
                zzlevbar(ig,1)=zphisbar(ig)/g
@@ -765 +762 @@
             endif         
          ENDDO
+
+      else !  if not moyzon
+      
+         DO l=1,nlayer   
+           zzlay_eff(ig,l)=g*rad*rad/(g*rad-(pphi(ig,l)+phisfi(ig)))-rad ! reference = geoid
+         ENDDO
+         zzlev_eff(ig,1)=phisfi(ig)/g
+         DO l=2,nlayer
+            z1=(pplay(ig,l-1)+pplev(ig,l))/ (pplay(ig,l-1)-pplev(ig,l))
+            z2=(pplev(ig,l)  +pplay(ig,l))/(pplev(ig,l)  -pplay(ig,l))
+            zzlev_eff(ig,l)=(z1*zzlay_eff(ig,l-1)+z2*zzlay_eff(ig,l))/(z1+z2)
+         ENDDO
+         zzlev_eff(ig,nlayer+1)=zzlay_eff(ig,nlayer)+(zzlay_eff(ig,nlayer)-zzlev_eff(ig,nlayer))
 
       endif  ! moyzon
@@ -1075 +1085 @@
             tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here
 #else
-            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
+            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi) ! JVO 19 : To be fixed, what altitude do we need ?
             pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
 #endif
@@ -1173 +1183 @@
 
                  ! Here we send zonal average fields ( corrected with cond ) from dynamics to chem. module
-                 call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,  &
+                 call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,zphisbar,  &
                               zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
                else ! 3D chemistry (or 1D run)
-                 call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,  &
-                              pplay,pplev,zzlay,zzlev,dycchi)
+                 call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,phisfi,  &
+                              pplay,pplev,zzlay_eff,zzlev_eff,dycchi)
                endif ! if moyzon
 
@@ -1602 +1612 @@
       CALL send_xios_field("dtrad",dtrad)
       CALL send_xios_field("dtdyn",zdtdyn)
+      CALL send_xios_field("dtdif",zdtdif)
 
       ! Chemical tracers