Changeset 3325


Ignore:
Timestamp:
May 13, 2024, 11:13:14 AM (8 months ago)
Author:
llange
Message:

Mars PCM
Some modifications in vdifc and pbl_parameters to include the effect of water buoyoncy on the sublimation of water ice.
Note: it only works with paleoclimate = .true. (since the model is not tuned with that ...).
LL

Location:
trunk/LMDZ.MARS
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.MARS/changelog.txt

    r3316 r3325  
    46264626== 25/04/2024 == JBC
    46274627Reversion of r3305 and r3307.
     4628
     4629== 13/05/2024 == LL
     4630Some modifications in vdifc and pbl_parameters to include the effect of water buoyoncy on the sublimation of water ice.
     4631Note: it only works with paleoclimate = .true. (since the model is not tuned with that ...).
  • trunk/LMDZ.MARS/deftank/field_def_physics_mars.xml

    r3305 r3325  
    542542             <field id="zdqsdif_ssi_tot"
    543543                   long_name="Total flux echanged with subsurface ice"
    544                    unit="kg.m-2.s-1" />       
     544                   unit="kg.m-2.s-1" />   
     545
     546
     547            <!-- Water ice sublimation parameters -->
     548            <field id="rib_dry_vdif_cd"
     549                   long_name="Dry Richardson number in vdif_cd"
     550                   unit="1" />
     551            <field id="rib_vdif_cd"
     552                   long_name="Richardson number in vdif_cd"
     553                   unit="1" />
     554            <field id="fm_vdif_cd"
     555                   long_name="Momentum stability function  in vdif_cd"
     556                   unit="1" />
     557            <field id="fh_vdif_cd"
     558                   long_name="Heat stability function  in vdif_cd"
     559                   unit="1" />
     560            <field id="z0t_vdif_cd"
     561                   long_name="Thermal roughness length  in vdif_cd"
     562                   unit="1" />   
     563            <field id="z0_vdif_cd"
     564                   long_name="Momentum roughness length  in vdif_cd"
     565                   unit="1" />   
     566            <field id="Reynolds_vdif_cd"
     567                   long_name="Reynolds number in vdif_cd"
     568                   unit="1" />   
    545569
    546570            <!-- CO2 cycle -->
  • trunk/LMDZ.MARS/libf/phymars/conf_phys.F

    r3272 r3325  
    3838      USE mod_phys_lmdz_transfert_para, ONLY: bcast
    3939      USE paleoclimate_mod,ONLY: paleoclimate,albedo_perennialco2,
    40      &                           lag_layer
     40     &                           lag_layer, include_waterbuoyancy
    4141      use microphys_h, only: mteta
    4242      use comsoil_h, only: adsorption_soil, choice_ads,ads_const_D,
     
    364364         lag_layer=.false.
    365365         call getin_p("lag_layer",lag_layer)
    366          write(*,*) " lag_layer = ", lag_layer
    367 
    368         write(*,*)"Is it paleoclimate run?"
     366         write(*,*) "lag_layer = ", lag_layer
     367
     368         write(*,*)"Is it paleoclimate run?"
    369369         paleoclimate=.false. ! default value
    370370         call getin_p("paleoclimate",paleoclimate)
    371          write(*,*)" paleoclimate = ",paleoclimate
    372 
     371         write(*,*)"paleoclimate = ",paleoclimate
    373372
    374373         write(*,*)"Albedo for perennial CO2 ice?"
     
    376375         call getin_p("albedo_perennialco2",albedo_perennialco2)
    377376         write(*,*)"albedo_perennialco2 = ",albedo_perennialco2
     377
     378         write(*,*)"Using lag layer??"
     379         lag_layer=.false.
     380         call getin_p("include_waterbuoyancy",include_waterbuoyancy)
     381         write(*,*) "include_waterbuoyancy = ", include_waterbuoyancy
     382
     383         if (include_waterbuoyancy.and.(.not.paleoclimate)) then
     384          print*,'You should not include the water buoyancy effect
     385     &        on current climate (model not tunned with this effect'
     386          call abort_physic(modname,
     387     &     "include_waterbuoyancy must be used with paleoclimate",1)
     388         endif
    378389
    379390! TRACERS:
  • trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90

    r3200 r3325  
    1515
    1616!$OMP THREADPRIVATE(paleoclimate)
    17     real,    save, allocatable, dimension(:,:) :: h2o_ice_depth       ! Thickness of the lag before H2O ice [m]
    18     real,    save, allocatable, dimension(:,:) :: lag_co2_ice         ! Thickness of the lag before CO2 ice [m]
    19     real,    save, allocatable, dimension(:,:) :: d_coef              ! Diffusion coeficent
    20     real,    save                              :: albedo_perennialco2 ! Albedo for perennial co2 ice [1]
    21     logical, save                              :: lag_layer           ! Does lag layer is present?
    22 !$OMP THREADPRIVATE(h2o_ice_depth,d_coef,lag_co2_ice,albedo_perennialco2)
     17    real,    save, allocatable, dimension(:,:) :: h2o_ice_depth           ! Thickness of the lag before H2O ice [m]
     18    real,    save, allocatable, dimension(:,:) :: lag_co2_ice             ! Thickness of the lag before CO2 ice [m]
     19    real,    save, allocatable, dimension(:,:) :: d_coef                  ! Diffusion coeficent
     20    real,    save                              :: albedo_perennialco2     ! Albedo for perennial co2 ice [1]
     21    logical, save                              :: lag_layer               ! Does lag layer is present?
     22    logical, save                              :: include_waterbuoyancy    ! Include the effect of water buoyancy when computing the sublimation of water ice ?
     23!$OMP THREADPRIVATE(h2o_ice_depth,d_coef,lag_co2_ice,albedo_perennialco2,include_waterbuoyancy)
    2324
    2425!=======================================================================
  • trunk/LMDZ.MARS/libf/phymars/pbl_parameters_mod.F90

    r3219 r3325  
    1919CONTAINS
    2020
    21       SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0,pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,tke,pts,ph,z_out,n_out, &
     21      SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0,pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,tke,pts,ph,pqvap,pqsurf,mumean,z_out,n_out, &
    2222                                T_out,u_out,ustar,tstar,vhf,vvv)
    2323                               
     
    2626      use turb_mod, only: turb_resolved
    2727      use lmdz_atke_turbulence_ini, only : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi
     28      use watersat_mod, only: watersat
     29      use paleoclimate_mod, only: include_waterbuoyancy
    2830
    2931      IMPLICIT NONE
     
    9294      REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)                 ! surface temperature, potentiel temperature
    9395      REAL, INTENT(IN) :: z_out(n_out)                              ! altitudes of the interpolation in the surface layer
     96      REAL, INTENT(IN) :: mumean(ngrid)                             ! Molecular mass of the atmosphere (kg/mol)
     97      REAL, INTENT(IN) :: pqvap(ngrid,nlay)                         ! Water vapor mixing ratio (kg/kg) to account for vapor flottability
     98      REAL, INTENT(IN) :: pqsurf(ngrid)                             ! Water ice frost on the surface (kg/m^2) to account for vapor flottability
    9499
    95100!    Outputs:
     
    142147      REAL f_ri_cd_min                                              ! minimum of the stability function with the ATKE scheme
    143148      REAL ric_4interp                                              ! critical richardson number which is either the one from Colaitis et al. (2013) or ATKE (1)
     149
     150
     151      REAL tsurf_v(ngrid)                                           ! Virtual surface temperature, accounting for vapor flottability
     152      REAL temp_v(ngrid)                                            ! Potential virtual air temperature in the frist layer, accounting for vapor flottability
     153      REAL mu_h2o                                                   ! Molecular mass of water vapor (kg/mol)
     154      REAL tol_frost                                                ! Tolerance to consider the effect of frost on the surface
     155      REAL qsat(ngrid)                                              ! saturated mixing ratio of water vapor
     156
     157
    144158!    Code:
    145159!    --------
     
    166180      tol_iter =  0.01
    167181      f_ri_cd_min = 0.01
     182      mu_h2o = 18e-3
     183      tol_frost = 1e-4
     184      tsurf_v(:) = 0.
     185      temp_v(:) = 0.
    168186
    169187! this formulation assumes alphah=1., implying betah=betam
     
    183201      endif
    184202     
    185      
     203
     204
     205      IF(include_waterbuoyancy) then
     206         temp_v(:) = ph(:,1)*(1.+pqvap(:,1)/(mu_h2o/mumean(:)))/(1.+pqvap(:,1))
     207         DO ig = 1,ngrid
     208            IF(pqsurf(ig).gt.tol_frost) then
     209               call watersat(1,pts(ig),pplay(ig,1),qsat(ig))
     210               tsurf_v(ig) = pts(ig)*(1.+qsat(ig)/(mu_h2o/mumean(ig)))/(1.+qsat(ig))
     211            ELSE
     212               tsurf_v(ig) = pts(ig)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1))
     213            ENDIF
     214         ENDDO
     215      ELSE
     216         tsurf_v(:) = pts(:)
     217         temp_v(:) =  ph(:,1)
     218      ENDIF
     219     
    186220      DO ig=1,ngrid
    187221         ite = 0.
     
    199233            IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then             
    200234            ! Richardson number formulation proposed by D.E. England et al. (1995)
    201                rib(ig) = (pg/pts(ig))*sqrt(zzlay(ig,1)*pz0(ig))*(((log(zzlay(ig,1)/pz0(ig)))**2)/(log(zzlay(ig,1)/pz0t)))*(ph(ig,1)-pts(ig))/zu2(ig)
     235               rib(ig) = (pg/tsurf_v(ig))*sqrt(zzlay(ig,1)*pz0(ig))*(((log(zzlay(ig,1)/pz0(ig)))**2)/(log(zzlay(ig,1)/pz0t)))*(temp_v(ig)-tsurf_v(ig))/zu2(ig)
    202236            ELSE
    203237               print*,'warning, infinite Richardson at surface'
  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r3304 r3325  
    25842584           call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
    25852585     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,q2,tsurf(:,iflat),
    2586      & zh,z_out,n_out,T_out,u_out,ustar,tstar,vhf,vvv)
     2586     & zh,zq(:,1,igcm_h2o_vap),qsurf(:,igcm_h2o_ice,iflat),mmean(:,1),
     2587     & z_out,n_out,T_out,u_out,ustar,tstar,vhf,vvv)
    25872588                   ! pourquoi ustar recalcule ici? fait dans vdifc.
    25882589
  • trunk/LMDZ.MARS/libf/phymars/vdif_cd_mod.F90

    r3219 r3325  
    1616CONTAINS
    1717
    18    SUBROUTINE vdif_cd(ngrid,nlay,nslope,pz0,pg,pz,pp,pu,pv,wstar,pts,ph,virtual,mumean,pqvap,pqsurf,pcdv,pcdh)
     18   SUBROUTINE vdif_cd(ngrid,nlay,nslope,pz0,pg,pz,pp,pu,pv,wstar,pts,ph,mumean,pqvap,pqsurf,write_outputs,pcdv,pcdh)
    1919   
    2020   use comcstfi_h
     
    2222   use watersat_mod, only: watersat
    2323   use lmdz_atke_turbulence_ini, only : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi
    24    
     24   use paleoclimate_mod, only: include_waterbuoyancy
     25   use write_output_mod, only: write_output
     26   use comslope_mod, ONLY: iflat
    2527   IMPLICIT NONE
    2628   include "callkeys.h"     
     29
     30
    2731!=======================================================================
    2832!
     
    4751!     pts(ngrid)          surface temperature
    4852!     ph(ngrid)           potential temperature T*(p/ps)^kappa
    49 !     virtual             Boolean to account for vapor flottability
    5053!     mumean              Molecular mass of the atmosphere (kg/mol)
    5154!     pqvap(ngrid,nlay)   Water vapor mixing ratio (kg/kg) to account for vapor flottability
     
    7679      REAL, INTENT(IN) :: pts(ngrid,nslope),ph(ngrid,nlay)  ! Surface Temperature and atmospheric temperature (K)
    7780      REAL, INTENT(IN) :: wstar(ngrid)                      ! Vertical velocity due to turbulence (m/s)
    78       LOGICAL, INTENT(IN) :: virtual                        ! Boolean to account for vapor flottability
    7981      REAL, INTENT(IN) :: mumean(ngrid)                     ! Molecular mass of the atmosphere (kg/mol)
    8082      REAL, INTENT(IN) :: pqvap(ngrid,nlay)                 ! Water vapor mixing ratio (kg/kg) to account for vapor flottability
    8183      REAL, INTENT(IN) :: pqsurf(ngrid,nslope)              ! Water ice frost on the surface (kg/m^2) to account for vapor flottability
     84      LOGICAL, INTENT(IN) :: write_outputs                  ! write_output in xios or not.
    8285     
    8386!   Outputs:
     
    100103!$OMP THREADPRIVATE(karman,nu)
    101104
    102       REAL rib(ngrid)            ! Bulk Richardson number (1)
    103       REAL fm(ngrid)             ! stability function for momentum (1)
    104       REAL fh(ngrid)             ! stability function for heat (1)
     105      REAL rib(ngrid,nslope)            ! Bulk Richardson number [virtual or dry] (1)
     106      REAL rib_dry(ngrid,nslope)        ! Dry bulk Richardson number [virtual or dry] (1)
     107      REAL fm(ngrid,nslope)             ! stability function for momentum (1)
     108      REAL fh(ngrid,nslope)             ! stability function for heat (1)
    105109
    106110! Formalism for stability functions  fm= 1/phim^2; fh = 1/(phim*phih) where
     
    111115      REAL betam, betah, alphah, bm, bh, lambda
    112116
    113       REAL pz0t                  ! initial thermal roughness length z0t for the iterative scheme (m)
    114       REAL ric_colaitis          ! critical richardson number (1)
    115       REAL reynolds(ngrid)       ! Reynolds number (1)
    116       REAL prandtl(ngrid)        ! Prandtl number  (1)
    117       REAL pz0tcomp(ngrid)       ! computed z0t during the iterative scheme (m)
    118       REAL ite                   ! Number of iteration (1)
    119       REAL itemax                ! Maximum number of iteration (1)
    120       REAL residual              ! Residual between pz0t and pz0tcomp (m)
    121       REAL tol_iter              ! Tolerance for the residual to ensure convergence (1=
    122 
    123       REAL zu2(ngrid)            ! Near surface winds (m/s)
    124      
    125       REAL cdn(ngrid)            ! neutral momentum  drag coefficient (1)
    126       REAL chn(ngrid)            ! neutral heat  drag coefficient (1)
    127       REAL z1z0,z1z0t            ! ratios z1/z0 and z1/z0T
    128       REAL z1,zcd0               ! Neutral roughness (m) and Cd/Ch coefficient when call richls is not called
    129       REAL tsurf_v(ngrid,nslope) ! Virtual surface temperature, accounting for vapor flottability
    130       REAL temp_v(ngrid)         ! Potential virtual air temperature in the frist layer, accounting for vapor flottability
    131       REAL mu_h2o                ! Molecular mass of water vapor (kg/mol)
    132       REAL tol_frost             ! Tolerance to consider the effect of frost on the surface
    133       REAL qsat(ngrid)           ! saturated mixing ratio of water vapor
     117      REAL pz0t                   ! initial thermal roughness length z0t for the iterative scheme (m)
     118      REAL ric_colaitis           ! critical richardson number (1)
     119      REAL reynolds(ngrid,nslope) ! Reynolds number (1)
     120      REAL prandtl(ngrid)         ! Prandtl number  (1)
     121      REAL pz0tcomp(ngrid)        ! computed z0t during the iterative scheme (m)
     122      REAL z0t(ngrid,nslope)      ! computed z0t at the last step the iterative scheme (m)
     123      REAL ite                    ! Number of iteration (1)
     124      REAL itemax                 ! Maximum number of iteration (1)
     125      REAL residual               ! Residual between pz0t and pz0tcomp (m)
     126      REAL tol_iter               ! Tolerance for the residual to ensure convergence (1=
     127
     128      REAL zu2(ngrid)             ! Near surface winds (m/s)
     129     
     130      REAL cdn(ngrid)             ! neutral momentum  drag coefficient (1)
     131      REAL chn(ngrid)             ! neutral heat  drag coefficient (1)
     132      REAL z1z0,z1z0t             ! ratios z1/z0 and z1/z0T
     133      REAL z1,zcd0                ! Neutral roughness (m) and Cd/Ch coefficient when call richls is not called
     134      REAL tsurf_v(ngrid,nslope)  ! Virtual surface temperature, accounting for vapor flottability
     135      REAL temp_v(ngrid)          ! Potential virtual air temperature in the frist layer, accounting for vapor flottability
     136      REAL mu_h2o                 ! Molecular mass of water vapor (kg/mol)
     137      REAL tol_frost              ! Tolerance to consider the effect of frost on the surface
     138      REAL qsat(ngrid)            ! saturated mixing ratio of water vapor
    134139           
    135       REAL sm                    ! Stability function computed with the ATKE scheme
    136       REAL f_ri_cd_min           ! minimum of the stability function with the ATKE scheme
     140      REAL sm                     ! Stability function computed with the ATKE scheme
     141      REAL f_ri_cd_min            ! minimum of the stability function with the ATKE scheme
    137142     
    138143!    Code:
     
    144149      mu_h2o = 18e-3
    145150      tol_frost = 1e-4
    146       reynolds(:) = 0.
     151      reynolds(:,:) = 0.
    147152      pz0t = 0.
    148153      pz0tcomp(:) = 0.1*pz0(:)
    149       rib(:) = 0.
     154      rib(:,:) = 0.
     155      rib_dry(:,:) = 0.
    150156      pcdv(:,:) = 0.
    151157      pcdh(:,:) = 0.
     158      z0t(:,:) = 0.
    152159      f_ri_cd_min = 0.01
    153160! this formulation assumes alphah=1., implying betah=betam
     
    162169      ric_colaitis = betah/(betam**2)
    163170     
    164       IF(virtual) then
     171      IF(include_waterbuoyancy) then
    165172         temp_v(:) = ph(:,1)*(1.+pqvap(:,1)/(mu_h2o/mumean(:)))/(1.+pqvap(:,1))
    166173         DO islope = 1,nslope
     
    210217                     IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.)
    211218! Richardson number formulation proposed by D.E. England et al. (1995)
    212                      rib(ig) = (pg/tsurf_v(ig,islope))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(temp_v(ig)-tsurf_v(ig,islope))/zu2(ig)
     219                     rib(ig,islope) = (pg/tsurf_v(ig,islope))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(temp_v(ig)-tsurf_v(ig,islope))/zu2(ig)
     220                     rib_dry(ig,islope) = (pg/pts(ig,islope))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(ph(ig,1)-pts(ig,islope))/zu2(ig)
    213221                   ELSE
    214222                      print*,'warning, infinite Richardson at surface'
    215223                      print*,pu(ig,1),pv(ig,1)
    216                       rib(ig) = ric_colaitis         
     224                      rib(ig,islope) = ric_colaitis   
     225                      rib_dry(ig,islope) = ric_colaitis
    217226                   ENDIF
    218227       
     
    221230                   IF(callatke) THEN
    222231                   ! Computation following parameterizaiton from ATKE
    223                       IF(rib(ig) .gt. 0) THEN
    224                          ! STABLE BOUNDARY LAYER :
    225                          sm = MAX(smmin,cn*(1.-rib(ig)/ric))
     232                      IF(rib(ig,islope) .gt. 0) THEN
     233                         ! STABLE BOUNDARY LAYER :include_waterbuoyancy
     234                         sm = MAX(smmin,cn*(1.-rib(ig,islope)/ric))
    226235                         ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
    227                          prandtl(ig) = pr_neut*exp(-pr_slope/pr_neut*rib(ig)+rib(ig)/pr_neut) + rib(ig) * pr_slope
    228                          fm(ig) = max((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
    229                          fh(ig) = max((fm(ig) / prandtl(ig)), f_ri_cd_min)
     236                         prandtl(ig) = pr_neut*exp(-pr_slope/pr_neut*rib(ig,islope)+rib(ig,islope)/pr_neut) + rib(ig,islope) * pr_slope
     237                         fm(ig,islope) = max((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig,islope) / prandtl(ig))**(1./2.)),f_ri_cd_min)
     238                         fh(ig,islope) = max((fm(ig,islope) / prandtl(ig)), f_ri_cd_min)
    230239                      ELSE
    231240                         ! UNSTABLE BOUNDARY LAYER :
    232                          sm = 2./rpi * (cinf - cn) * atan(-rib(ig)/ri0) + cn
    233                          prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig)/ri1) + pr_neut
    234                          fm(ig) = MAX((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
    235                          fh(ig) = MAX((fm(ig) / prandtl(ig)), f_ri_cd_min)
     241                         sm = 2./rpi * (cinf - cn) * atan(-rib(ig,islope)/ri0) + cn
     242                         prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig,islope)/ri1) + pr_neut
     243                         fm(ig,islope) = MAX((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig,islope) / prandtl(ig))**(1./2.)),f_ri_cd_min)
     244                         fh(ig,islope) = MAX((fm(ig,islope) / prandtl(ig)), f_ri_cd_min)
    236245                      ENDIF ! Rib < 0
    237246                   ELSE
    238247                   ! Computation following parameterizaiton from from D.E. England et al. (95)
    239                       IF (rib(ig) .gt. 0.) THEN
     248                      IF (rib(ig,islope) .gt. 0.) THEN
    240249                        ! STABLE BOUNDARY LAYER :
    241250                        prandtl(ig) = 1.
    242                         IF(rib(ig) .lt. ric_colaitis) THEN
     251                        IF(rib(ig,islope) .lt. ric_colaitis) THEN
    243252               ! Assuming alphah=1. and bh=bm for stable conditions :
    244                            fm(ig) = ((ric_colaitis-rib(ig))/ric_colaitis)**2
    245                            fh(ig) = ((ric_colaitis-rib(ig))/ric_colaitis)**2
     253                           fm(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2
     254                           fh(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2
    246255                        ELSE
    247256               ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
    248                            fm(ig) = 1.
    249                            fh(ig) = 1.
     257                           fm(ig,islope) = 1.
     258                           fh(ig,islope) = 1.
    250259                         ENDIF
    251260                     ELSE
    252261                        ! UNSTABLE BOUNDARY LAYER :
    253                         fm(ig) = sqrt(1.-lambda*bm*rib(ig))
    254                         fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25
    255                         prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5)
     262                        fm(ig,islope) = sqrt(1.-lambda*bm*rib(ig,islope))
     263                        fh(ig,islope) = (1./alphah)*((1.-lambda*bh*rib(ig,islope))**0.5)*(1.-lambda*bm*rib(ig,islope))**0.25
     264                        prandtl(ig) = alphah*((1.-lambda*bm*rib(ig,islope))**0.25)/((1.-lambda*bh*rib(ig,islope))**0.5)
    256265                     ENDIF ! Rib < 0
    257266                  ENDIF ! atke
    258267              ! Recompute the reynolds and z0t
    259                   reynolds(ig) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu)
    260                   pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig)**0.25)*(prandtl(ig)**0.5)+5*karman)
     268                  reynolds(ig,islope) = karman*sqrt(fm(ig,islope))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu)
     269                  pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig,islope)**0.25)*(prandtl(ig)**0.5)+5*karman)
    261270                  residual = abs(pz0t-pz0tcomp(ig))
    262271                  ite = ite+1
     
    264273           
    265274! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions       
    266             pcdv(ig,islope)=cdn(ig)*fm(ig)
    267             pcdh(ig,islope)=chn(ig)*fh(ig)         
     275            pcdv(ig,islope)=cdn(ig)*fm(ig,islope)
     276            pcdh(ig,islope)=chn(ig)*fh(ig,islope)
     277            z0t(ig,islope) =  pz0tcomp(ig)       
    268278            pz0t = 0. ! for next grid point
    269279            ENDDO ! of ngrid
    270280         enddo
    271281      ENDIF !of if call richardson surface layer
     282     
     283      IF(write_outputs) then
     284
     285         call write_output("rib_dry_vdif_cd","Dry Richardson number in vdif_cd","1",rib_dry(:,iflat))
     286
     287         call write_output("rib_vdif_cd","Richardson number in vdif_cd","1",rib(:,iflat))
     288
     289         call write_output("fm_vdif_cd","Momentum stability function  in vdif_cd","1",fm(:,iflat))
     290
     291         call write_output("fh_vdif_cd","Heat stability function  in vdif_cd","1",fh(:,iflat))
     292
     293         call write_output("z0t_vdif_cd","Thermal roughness length  in vdif_cd","m",z0t(:,iflat))
     294
     295         call write_output("z0_vdif_cd","Momentum roughness length  in vdif_cd","m",pz0(:))
     296
     297         call write_output("Reynolds_vdif_cd","Reynolds number in vdif_cd","1",reynolds(:,iflat))
     298
     299      ENDIF
    272300
    273301      RETURN
  • trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F

    r3295 r3325  
    235235      LOGICAL :: writeoutput          ! boolean to say to soilexchange.F if we are at the last iteration and thus if he  can write in the diagsoil
    236236      REAL, INTENT(out) :: pore_icefraction(ngrid,nsoil,nslope) ! ice filling fraction in the pores
    237 !! Water buyoncy
    238       LOGICAL :: virtual
    239237
    240238!! Subsurface ice interactions
     
    309307      zq1temp_regolith(1:ngrid)=0
    310308      zdqsdif_tot(1:ngrid)=0
    311       virtual = .false.
    312309      pore_icefraction(:,:,:) = 0.
    313310      zdqsdif_ssi_atm(:,:) = 0.
     
    434431
    435432      CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,pu,pv,wstar,
    436      &          ptsrf,ph,virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),
    437      &          pqsurf(:,igcm_h2o_ice,:),
     433     &          ptsrf,ph,mmean(:,1),zq(:,:,igcm_h2o_vap),
     434     &          pqsurf(:,igcm_h2o_ice,:), .false.,
    438435     &          zcdv_true,zcdh_true)
    439436
     
    612609
    613610      CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,zu,zv,wstar,
    614      &          ptsrf,ph,virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),
    615      &          pqsurf(:,igcm_h2o_ice,:),
     611     &          ptsrf,ph,mmean(:,1),zq(:,:,igcm_h2o_vap),
     612     &          pqsurf(:,igcm_h2o_ice,:), .true.,
    616613     &          zcdv_true,zcdh_true)
    617614
Note: See TracChangeset for help on using the changeset viewer.