Changeset 3153 for trunk/LMDZ.MARS

Timestamp:

Dec 9, 2023, 6:13:21 PM (12 months ago)

Author:

llange

Message:

Mars PCM
Add the possibility to compute Cd; Ch based on the virtual potential temperature to account for water flotability.
To do so, a boolean "virtual", set to false by default, must be set to true (for now hard coded as future modifications will follow).

LL

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/vdif_cd_mod.F90 (modified) (7 diffs)
• libf/phymars/vdifc_mod.F (modified) (5 diffs)

trunk/LMDZ.MARS/changelog.txt

@@ -4388 +4388 @@
 
 == 07/12/2023 == LL
-
 Switching pbl_parameter.F into a module. Cleaning of unused variables/sections in it.
 Correction in the computations of ustar, thetastar as it was done in a different way compare to vdif_cd. It is now self-consistent.
 Maybe pbl_surface in the MCD needs also to be changed ? In discussion with Ehouarn and Aymeric.
+
+== 09/12/2023 == LL
+Add the possibility to compute Cd; Ch based on the virtual potential temperature to account for water flotability.
+To do so, a boolean "virtual", set to false by default, must be set to true (for now hard coded as future modifications will follow).

trunk/LMDZ.MARS/libf/phymars/vdif_cd_mod.F90

@@ -16 +16 @@
 CONTAINS
 
-   SUBROUTINE vdif_cd(ngrid,nlay,pz0,pg,pz,pu,pv,wstar,pts,ph,pcdv,pcdh)
+   SUBROUTINE vdif_cd(ngrid,nlay,pz0,pg,pz,pp,pu,pv,wstar,pts,ph,virtual,mumean,pqvap,pqsurf,pcdv,pcdh)
    
    use comcstfi_h
    use turb_mod, only: turb_resolved
-
+   use watersat_mod, only: watersat
 
    IMPLICIT NONE
@@ -38 +38 @@
 !   inputs:
 !   ------
-!     ngrid            size of the horizontal grid
-!     pg               gravity (m s -2)
-!     pz(ngrid,nlay)   height of layers
-!     pu(ngrid,nlay)   u component of the wind
-!     pv(ngrid,nlay)   v component of the wind
-!     pts(ngrid)       surface temperature
-!     ph(ngrid)        potential temperature T*(p/ps)^kappa
+!     ngrid               size of the horizontal grid
+!     pg                  gravity (m s -2)
+!     pz(ngrid,nlay)      height of layers
+!     pp(ngrid,nlay)      pressure of layers
+!     pu(ngrid,nlay)      u component of the wind
+!     pv(ngrid,nlay)      v component of the wind
+!     pts(ngrid)          surface temperature
+!     ph(ngrid)           potential temperature T*(p/ps)^kappa
+!     virtual             Boolean to account for vapor flottability
+!     mumean              Molecular mass of the atmosphere (kg/mol)
+!     pqvap(ngrid,nlay)   Water vapor mixing ratio (kg/kg) to account for vapor flottability
+!     pqsurf(ngrid)       Water ice frost on the surface(kg/m^2) to account for vapor flottability
 !
 !   outputs:
@@ -67 +72 @@
       INTEGER, INTENT(IN) :: ngrid,nlay                     ! Number of points in the physical grid and atmospheric grid
       REAL, INTENT(IN) :: pz0(ngrid)                        ! Surface Roughness (m)
-      REAL, INTENT(IN) :: pg,pz(ngrid,nlay)                 ! Mars gravity (m/s^2)
+      REAL, INTENT(IN) :: pg                                ! Mars gravity (m/s^2)
+      REAL, INTENT(IN) :: pz(ngrid,nlay),pp(ngrid,nlay)     ! Layers pseudo altitudes (m) and pressure (Pa)                
       REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)     ! Zonal and Meriditionnal  winds (m/s)
       REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)         ! Surface Temperature and atmospheric temperature (K)
       REAL, INTENT(IN) :: wstar(ngrid)                      ! Vertical velocity due to turbulence (m/s)
+      LOGICAL, INTENT(IN) :: virtual                        ! Boolean to account for vapor flottability
+      REAL, INTENT(IN) :: mumean(ngrid)                     ! Molecular mass of the atmosphere (kg/mol)
+      REAL, INTENT(IN) :: pqvap(ngrid,nlay)                 ! Water vapor mixing ratio (kg/kg) to account for vapor flottability
+      REAL, INTENT(IN) :: pqsurf(ngrid)                     ! Water ice frost on the surface (kg/m^2) to account for vapor flottability
       
 !   Outputs:
@@ -117 +127 @@
       REAL z1z0,z1z0t        ! ratios z1/z0 and z1/z0T
       REAL z1,zcd0           ! Neutral roughness (m) and Cd/Ch coefficient when call richls is not called 
-      
+      REAL tsurf_v(ngrid)    ! Virtual surface temperature, accounting for vapor flottability
+      REAL temp_v(ngrid)     ! Potential virtual air temperature in the frist layer, accounting for vapor flottability
+      REAL mu_h2o            ! Molecular mass of water vapor (kg/mol)
+      REAL tol_frost         ! Tolerance to consider the effect of frost on the surface
+      REAL qsat(ngrid)       ! saturated mixing ratio of water vapor 
+           
 !    Code:
 !    --------
@@ -124 +139 @@
       itemax= 10 
       tol_iter =  0.01
+      mu_h2o = 18e-3
+      
       reynolds(:) = 0.
       pz0t = 0.
@@ -141 +158 @@
       
       ric = betah/(betam**2)
+      
+      IF(virtual) then
+         DO ig = 1,ngrid
+            temp_v(ig) = ph(ig,1)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1))
+            call watersat(1,pts(ig),pp(ig,1),qsat(ig))
+            tsurf_v(ig) = pts(ig)*(1.+qsat(ig)/(mu_h2o/mumean(ig)))/(1.+qsat(ig))
+         ENDDO
+      ELSE
+         tsurf_v(:) = pts(:)
+         temp_v(:) =  ph(:,1)
+      ENDIF
 
       IF(.not.callrichsl) then
@@ -170 +198 @@
                   IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.) 
 ! Richardson number formulation proposed by D.E. England et al. (1995)
-                  rib(ig) = (pg/pts(ig))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(ph(ig,1)-pts(ig))/zu2(ig)
+                  rib(ig) = (pg/tsurf_v(ig))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(temp_v(ig)-tsurf_v(ig))/zu2(ig)
                 ELSE
                    print*,'warning, infinite Richardson at surface'

trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F

@@ -208 +208 @@
       REAL,SAVE :: m_co2, m_noco2, A , B
       REAL vmr_co2(ngrid,nlay)
-      REAL qco2,mmean
+      REAL qco2,mmean(ngrid,nlay)
 
 !$OMP THREADPRIVATE(ico2,m_co2,m_noco2,A,B)
@@ -235 +235 @@
       REAL zdqsdif_tot(ngrid)         ! subtimestep pdqsdif for water ice
       LOGICAL :: writeoutput          ! boolean to say to soilexchange.F if we are at the last iteration and thus if he  can write in the diagsoil
+      
+!! Water buyoncy
+      LOGICAL :: virtual
+
+
 
 c    ** un petit test de coherence
@@ -297 +302 @@
       zdqsdif_tot(1:ngrid)=0
       h2o_ice_depth(1:ngrid,1:nslope)=1
+      
+      virtual = .false.
+      
 c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
 c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
@@ -354 +362 @@
      &             ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep)
 c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
-              mmean=1/(A*qco2 +B)
-              vmr_co2(ig,ilev) = qco2*mmean/m_co2
+              mmean(ig,ilev)=1/(A*qco2 +B)
+              vmr_co2(ig,ilev) = qco2*mmean(ig,ilev)/m_co2
             ENDDO
          ENDDO
@@ -412 +420 @@
 c       ---------------------
 
-      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf_tmp
-     &          ,ph,zcdv_true,zcdh_true)
+      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pplay,pu,pv,wstar,ptsrf_tmp,
+     &          ph,virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),
+     &          pqsurf(:,igcm_h2o_ice,major_slope(:)), 
+     &          zcdv_true,zcdh_true)
 
         zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)