Changeset 4044 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Feb 5, 2026, 3:39:30 PM (3 weeks ago)
Author:
jmauxion
Message:

Mars PCM:
improvedclouds: from now on, the tendencies from previous physics are added
all-at-once rather than distributed over the microphysics substeps.
JM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F90

    r3939 r4044  
    7979REAL, dimension(ngrid,nlay) :: alpha ! HDO equilibrium fractionation coefficient (Saturation=1)
    8080REAL, dimension(ngrid,nlay) :: alpha_c ! HDO real fractionation coefficient
    81 ! REAL :: sumcheck
    8281REAL*8 :: ph2o ! Water vapor partial pressure (Pa)
    8382REAL*8 :: satu ! Water vapor saturation ratio over ice
    … …  
    9392REAL :: Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
    9493
    95 !     Parameters of the size discretization used by the microphysical scheme
     94! Parameters of the size discretization used by the microphysical scheme
    9695DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
    9796DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
    … …  
    111110REAL, dimension(ngrid,nlay) :: subpdtcloud ! Temperature variation due to latent heat
    112111REAL, dimension(ngrid,nlay,nq) :: zq0 ! local initial value of tracers
    113 ! REAL, dimension(ngrid,nlay,nq) :: zt0 ! local initial value of temperature
    114112INTEGER, dimension(ngrid,nlay) :: zimicro ! Subdivision of ptimestep
    115113INTEGER, dimension(ngrid,nlay) :: count_micro ! Number of microphys calls
    116 REAL, dimension(ngrid,nlay) :: zpotcond_inst ! condensable water at the beginning of physics (kg/kg)
    117 REAL, dimension(ngrid,nlay) :: zpotcond_full ! condensable water with integrated tendancies (kg/kg)
    118114REAL, dimension(ngrid,nlay) :: zpotcond ! maximal condensable water (previous two)
    119115REAL, dimension(1) :: zqsat_tmp ! maximal condensable water (previous two)
    … …  
    154150  rad_cld(1) = rmin_cld
    155151  vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
    156   ! vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld
    157152
    158153  do i=1,nbin_cld-1
    … …  
    179174
    180175  ! we save that so that it is not computed at each timestep and gridpoint
    181   ! rb_cld(i) = log(rb_cld(i)) 
    182176  rb_cld(:) = log(rb_cld(:))
    183177
    … …  
    207201rice(:,:) = 1.e-8
    208202
    209 !     Initialize the temperature, tracers
    210 zt(:,:)=pt(:,:)
    211 zq(:,:,:)=pq(:,:,:)
     203! Initialize the temperature and tracers with tendencies
     204
     205! If scavenging, add tendency for dust all-at-once
     206IF (scavenging) THEN
     207   zq(:,:,igcm_dust_mass) =zq(:,:,igcm_dust_mass)+pdq(:,:,igcm_dust_mass)*ptimestep
     208   zq(:,:,igcm_dust_number) =zq(:,:,igcm_dust_number)+pdq(:,:,igcm_dust_number)*ptimestep
     209ENDIF ! scavenging
     210
     211! Add tendency for ccn all-at-once
     212zq(:,:,igcm_ccn_mass) = zq(:,:,igcm_ccn_mass) + pdq(:,:,igcm_ccn_mass)*ptimestep     
     213zq(:,:,igcm_ccn_number) = zq(:,:,igcm_ccn_number) + pdq(:,:,igcm_ccn_number)*ptimestep
     214
     215! Add tendency for water all-at-once
     216zq(:,:,igcm_h2o_ice) = zq(:,:,igcm_h2o_ice)+pdq(:,:,igcm_h2o_ice)*ptimestep
     217zq(:,:,igcm_h2o_vap) = zq(:,:,igcm_h2o_vap)+pdq(:,:,igcm_h2o_vap)*ptimestep
     218
     219! Add tendency for HDO (if computed) all-at-once
     220IF (hdo) THEN
     221   zq(:,:,igcm_hdo_ice) = zq(:,:,igcm_hdo_ice)+pdq(:,:,igcm_hdo_ice)*ptimestep
     222   zq(:,:,igcm_hdo_vap) = zq(:,:,igcm_hdo_vap)+pdq(:,:,igcm_hdo_vap)*ptimestep
     223ENDIF
     224
     225! Add tendency for temp all-at-one
     226zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
     227
     228! Local temp tendency from clouds (due to latent heath release)
    212229subpdtcloud(:,:)=0
    213230
     231! Handle very small values to prevent precision error
    214232WHERE( zq(:,:,:) < 1.e-30 ) zq(:,:,:) = 1.e-30
    215233
    … …  
    221239
    222240! Compute the condensable amount of water vapor or the sublimable water ice (if negative)
    223 ! Attention ici pdt*ptimestep peut etre severe
    224 call watersat(ngrid*nlay,max(1.,zt+pdt*ptimestep),pplay,zqsat) ! DIFF: "temp+trend" at least 1
    225 zpotcond_full=(zq(:,:,igcm_h2o_vap)+pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat
    226 zpotcond_full=max(zpotcond_full,-zq(:,:,igcm_h2o_ice))
    227 call watersat(ngrid*nlay,zt,pplay,zqsat)
    228 zpotcond_inst=zq(:,:,igcm_h2o_vap) - zqsat
    229 zpotcond_inst=max(zpotcond_inst,-zq(:,:,igcm_h2o_ice))
    230 ! zpotcond is the most strict criterion between the two
    231 DO l=1,nlay
    232   DO ig=1,ngrid
    233     if (zpotcond_full(ig,l).gt.0.) then
    234       zpotcond(ig,l)=max(zpotcond_inst(ig,l),zpotcond_full(ig,l))
    235     else if (zpotcond_full(ig,l).le.0.) then
    236       zpotcond(ig,l)=min(zpotcond_inst(ig,l),zpotcond_full(ig,l))
    237     endif ! (zpotcond_full.gt.0.) then
    238     ! zpotcond(ig,l)=1.e-2
    239   ENDDO !ig=1,ngrid
    240 ENDDO !l=1,nlay
     241call watersat(ngrid*nlay,max(1.,zt),pplay,zqsat) ! Make sure "temp+tendency" at least 1
     242zpotcond=zq(:,:,igcm_h2o_vap) - zqsat
    241243     
    242244countcells = 0
    … …  
    250252    IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l), zimicro(ig,l))
    251253    spenttime = 0.
    252     dMicetot= 0. ! DIFF: dMicetot=new var
     254    dMicetot= 0.
    253255    ending_ts=.false.
    254256    DO while (.not.ending_ts)
    … …  
    259261      satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
    260262      microtimestep=ptimestep/real(zimicro(ig,l))
    261       ! if (satu.lt.1.0) then
    262       !   microtimestep=ptimestep/30.
    263       !   write(*,*) "saturation correction" !JN
    264       ! endif
    265263     
    266264      ! Initialize tracers for scavenging + hdo computations (JN)
    … …  
    298296        enddo
    299297     
    300         ! sumcheck = 0
    301         ! do i = 1, nbin_cld
    302         !   sumcheck = sumcheck + n_aer(i)
    303         ! enddo
    304         ! sumcheck = abs(sumcheck/No - 1)
    305         ! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
    306         !   print*, "WARNING, No sumcheck PROBLEM"
    307         !   print*, "sumcheck, No",sumcheck, No
    308         !   print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
    309         !   print*, "Dust binned distribution", n_aer
    310         ! endif
    311         !       
    312         ! sumcheck = 0
    313         ! do i = 1, nbin_cld
    314         !   sumcheck = sumcheck + m_aer(i)
    315         ! enddo
    316         ! sumcheck = abs(sumcheck/Mo - 1)
    317         ! if ((sumcheck .gt. 1e-5) .and.  (1./Rn .gt. rmin_cld)) then
    318         !   print*, "WARNING, Mo sumcheck PROBLEM"
    319         !   print*, "sumcheck, Mo",sumcheck, Mo
    320         !   print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
    321         !   print*, "Dust binned distribution", m_aer
    322         ! endif
    323 
    324298        ! Get the rates of nucleation
    325299        call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
    … …  
    368342
    369343        ! With the above scheme, dMice cannot be bigger than vapor, but can be bigger than all available ice.
    370         ! dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
    371         dMice = max(dMice,-zq(ig,l,igcm_h2o_ice) - pdq(ig,l,igcm_h2o_ice)*microtimestep) ! DIFF: take trend into account
    372         ! dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
    373         dMice = min(dMice,zq(ig,l,igcm_h2o_vap) + pdq(ig,l,igcm_h2o_vap)*microtimestep)
     344        dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
     345        dMice = min(dMice,zq(ig,l,igcm_h2o_vap))
    374346
    375347        zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
    … …  
    400372              ! Calculation of the fractionnation coefficient at equilibrium
    401373              alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
    402               ! alpha = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
    403374              ! Calculation of the 'real' fractionnation coefficient
    404375              alpha_c(ig,l) = (alpha(ig,l)*satu)/( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
    405               ! alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
    406376            else
    407377              alpha_c(ig,l) = 1.d0
    … …  
    459429      ! No more getting tendency : we increment tracers & temp directly
    460430
    461       ! Increment tracers & temp for following microtimestep
    462       ! Dust :
    463       ! Special treatment of dust_mass / number for scavenging ?
    464       IF (scavenging) THEN
    465         zq(ig,l,igcm_dust_mass) =zq(ig,l,igcm_dust_mass)+pdq(ig,l,igcm_dust_mass)*microtimestep
    466         zq(ig,l,igcm_dust_number) =zq(ig,l,igcm_dust_number)+pdq(ig,l,igcm_dust_number)*microtimestep
    467       ELSE
    468         zq(ig,l,igcm_dust_mass) =zq0(ig,l,igcm_dust_mass)
    469         zq(ig,l,igcm_dust_number) =zq0(ig,l,igcm_dust_number)
    470       ENDIF !(scavenging) THEN
    471       zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + pdq(ig,l,igcm_ccn_mass)*microtimestep
    472       zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)*microtimestep
    473 
    474       ! Water :
    475       zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*microtimestep
    476       zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*microtimestep
    477 
    478       ! HDO (if computed) :
    479       IF (hdo) THEN
    480         zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+pdq(ig,l,igcm_hdo_ice)*microtimestep
    481         zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)+pdq(ig,l,igcm_hdo_vap)*microtimestep
    482       ENDIF ! hdo
    483 
    484       ! Could also set subpdtcloud to 0 if not activice to make it simpler or change name of the flag
     431      ! If not activice, the tendency from latent heat release is set to zero
    485432      IF (.not.activice) subpdtcloud(ig,l)=0.
    486433
    487       ! Temperature :
    488       zt(ig,l) = zt(ig,l)+(pdt(ig,l)+subpdtcloud(ig,l))*microtimestep
     434      ! Temperature change as a feedback from latent heat release
     435      zt(ig,l) = zt(ig,l)+subpdtcloud(ig,l)*microtimestep
    489436
    490437      ! Prevent negative tracers ! JN
    … …  
    493440      IF (cloud_adapt_ts) then
    494441        ! Estimation of how much is actually condensing/subliming
    495          dMicetot=dMicetot+abs(dMice) ! DIFF: accumulation of absolute dMice
    496         ! dMicetot=dMicetot+dMice
    497         ! write(*,*) "dMicetot= ", dMicetot
    498         ! write(*,*) "we are in (l,ig):", l,ig !JN
     442         dMicetot=dMicetot+abs(dMice) ! total accumulated ice formation since the beginning
    499443        IF (spenttime.ne.0) then
    500444          zdq=(dMicetot/spenttime)!*(ptimestep-spenttime)
    501         ELSE
    502           !Initialization for spenttime=0
    503           zdq=zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)
    504445        ENDIF
    505         ! Threshold with powerlaw (sanity check)
    506         ! zdq=min(abs(zdq),
    507         ! & abs(zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)))
    508446        zdq=abs(zdq)
    509447        call adapt_imicro(ptimestep,zdq,zimicro(ig,l))
    … …  
    517455
    518456!------ Useful outputs to check how it went
    519 call write_output("zpotcond_inst","zpotcond_inst microphysics","(kg/kg)",zpotcond_inst(:,:))
    520 call write_output("zpotcond_full","zpotcond_full microphysics","(kg/kg)",zpotcond_full(:,:))
    521457call write_output("zpotcond","zpotcond microphysics","(kg/kg)",zpotcond(:,:))
    522458call write_output("count_micro","count_micro after microphysics","integer",count_micro(:,:))
    … …  
    524460!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
    525461IF (test_flag) then
    526   ! error2d(:) = 0.
    527462  DO l=1,nlay
    528463    DO ig=1,ngrid
    529       ! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
    530464      satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    531465      satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    532466    ENDDO
    533467  ENDDO
    534 
    535468  print*, 'count is ',countcells, ' i.e. ', countcells*100/(nlay*ngrid), '% for microphys computation'
    536469ENDIF ! endif test_flag
    537470
    538471END SUBROUTINE improvedclouds
    539 
    540 !=============================================================
    541 ! The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
    542 ! It is an analytical solution to the ice radius growth equation,
    543 ! with the approximation of a constant 'reduced' cunningham correction factor
    544 ! (lambda in growthrate.F) taken at radius req instead of rice   
    545 !=============================================================
    546 
    547 ! subroutine phi(rice,req,coeff1,coeff2,time)
    548 !     
    549 ! implicit none
    550 !     
    551 ! ! inputs
    552 ! real rice ! ice radius
    553 ! real req  ! ice radius at equilibirum
    554 ! real coeff1  ! coeff for the log
    555 ! real coeff2  ! coeff for the arctan
    556 !
    557 ! ! output     
    558 ! real time
    559 !     
    560 ! !local
    561 ! real var
    562 !     
    563 !  1.73205 is sqrt(3)
    564 !     
    565 ! var = max(
    566 ! &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
    567 !           
    568 ! time =
    569 ! &   coeff1 *
    570 ! &   log( var )
    571 ! & + coeff2 * 1.73205 *
    572 ! &   atan( (2*rice+req) / (1.73205*req) )
    573 !     
    574 ! return
    575 ! end
    576 
    577 !=======================================================================
    578472
    579473SUBROUTINE adapt_imicro(ptimestep,potcond,zimicro)
    … …  
    598492! alpha=1.81846504e+11
    599493! beta=1.54550140e+00
    600 alpha=1.88282793e+05 ! DIFF: new power law coefficients
    601 beta=4.57764370e-01
    602 zimicro=ceiling(coef*min(max(alpha*abs(potcond)**beta,50.),7000.))
    603 zimicro=2*zimicro ! DIFF: prefiction times two, allow to complete year
     494alpha=1.88282793e+05 ! Latest values for high obliquity
     495beta=4.57764370e-01  ! Latest values for high obliquity
     496!alpha=1.72198978e+10 ! Present day Mars
     497!beta=1.88473210e+00 
     498zimicro=ceiling(coef*min(max(alpha*abs(potcond)**beta,5.),7000.))
     499!zimicro=2*zimicro ! Prediction times two, allow to complete year at high obliquity
    604500
    605501END SUBROUTINE adapt_imicro
Note: See TracChangeset for help on using the changeset viewer.