Changeset 4724 for LMDZ6


Ignore:
Timestamp:
Oct 11, 2023, 9:27:08 AM (15 months ago)
Author:
evignon
Message:

debut du travail de replayisation des routines de neige soufflee
+ correction bug pour les temperatures positives

Location:
LMDZ6/trunk/libf
Files:
6 edited
6 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90

    r4723 r4724  
    1 module blowing_snow_ini_mod
     1module lmdz_blowing_snow_ini
    22
    33implicit none
    … …  
    1515   !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs)
    1616   !$OMP THREADPRIVATE(iflag_saltation_bs)
     17
     18    real tbsmelt                  ! parameter to calculate melting fraction of BS sedimentation
     19    parameter (tbsmelt=278.15)
     20
     21    real taumeltbs0               ! Melting time scale of blowing snow at 273.15K
     22    parameter (taumeltbs0=1800.0)
     23
     24    real qbsmin                   ! Minimum blowing snow specific content
     25    parameter (qbsmin=1.E-10)
     26
    1727
    1828      contains
    … …  
    7080      end subroutine blowing_snow_ini
    7181
    72 end module blowing_snow_ini_mod
     82end module lmdz_blowing_snow_ini

  • LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90

    r4723 r4724  
    1 subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
     1module lmdz_blowing_snow_sublim_sedim
     2
     3contains
     4subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
    25
    36!==============================================================================
    … …  
    811
    912
    10 use blowing_snow_ini_mod, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs
    11 use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2
     13use lmdz_blowing_snow_ini, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs, qbsmin
     14use lmdz_blowing_snow_ini, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, tbsmelt, taumeltbs0
    1215USE lmdz_lscp_tools, only : calc_qsat_ecmwf
    1316
    … …  
    2528real, intent(in)                        :: dtime
    2629real, intent(in), dimension(ngrid,nlay) :: temp
    27 real, intent(in), dimension(ngrid,nlay) :: q
     30real, intent(in), dimension(ngrid,nlay) :: qv
    2831real, intent(in), dimension(ngrid,nlay) :: qbs
    2932real, intent(in), dimension(ngrid,nlay) :: pplay
    … …  
    4952integer                                  :: k,i,n
    5053real                                     :: zqev0, zqevi, zqevti, zcpair, zcpeau, dqbsmelt
    51 real                                     :: dqsedim,precbs
     54real                                     :: dqsedim,precbs, deltaqchaud, zmelt, taumeltbs
     55real                                     :: maxdeltaqchaud
    5256real, dimension(ngrid)                   :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn
    5357real, dimension(ngrid)                   :: zqbsi,zmqc, zmair, zdz
    … …  
    7377    DO i=1,ngrid
    7478        zt(i)=temp(i,k)
    75         zq(i)=q(i,k)
     79        zq(i)=qv(i,k)
    7680        zqbs(i)=qbs(i,k)
    7781    ENDDO
    … …  
    111115   
    112116    DO i = 1, ngrid
    113        ! if sedimentation:
    114        IF (sedim(i) .GT. 0.) THEN
     117
     118       zrho(i,k)  = pplay(i,k) / zt(i) / RD
     119       zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
     120       ! BS fall velocity
     121       velo(i,k) = fallv_bs
     122
    115123
    116124          IF (zt(i) .GT. RTT) THEN
     125             
    117126             ! if positive celcius temperature, we assume
    118              ! that all the blowing snow melt and evaporate
    119               zqevti=sedim(i)*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG     
    120              ! we ensure that the whole mesh does not exceed saturation wrt liquid
    121               zqev0 = MAX(0.0, qsl(i)-zq(i))
    122               zqevi = MIN(zqev0,zqevti)
    123               ! New solid precipitation fluxes
    124               sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
    125               ! vapor, temperature, precip fluxes update
    126               zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    127               zq(i) = max(0., zq(i))
    128               !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    129               !zqbs(i) = max(0., zqbs(i))
    130               ! melting
    131               zt(i) = zt(i)                             &
    132                 + (sedimn(i)-sedim(i))                      &
    133                 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
    134                 * RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
    135               ! evaporation
    136               zt(i) = zt(i)                             &
    137                 + (sedimn(i)-sedim(i))                      &
    138                 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
    139                 * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))     
     127             ! that part of the the blowing snow flux melts and evaporates
     128             
     129             ! vapor, bs, temperature, precip fluxes update
     130             zmelt = ((zt(i)-RTT)/(tbsmelt-RTT))
     131             zmelt = MIN(MAX(zmelt,0.),1.)
     132             sedimn(i)=sedim(i)*zmelt
     133             deltaqchaud=-(sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     134             ! max evap such as celcius temperature cannot become negative
     135             maxdeltaqchaud= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
     136
     137             deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud)
     138             zq(i) = zq(i) + deltaqchaud 
     139             
     140             ! melting + evaporation
     141             zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
     142
     143             sedim(i)=sedimn(i)
     144
     145             ! if temperature still positive, we assume that part of the blowing snow
     146             ! already present in the mesh melts and evaporates with a typical time
     147             ! constant between taumeltbs0 and 0
     148 
     149             IF ( zt(i) .GT. RTT ) THEN
     150                taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT)))
     151                deltaqchaud=min(zqbs(i),zqbs(i)/taumeltbs*dtime)
     152                maxdeltaqchaud= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
     153                deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud)
     154                zq(i) = zq(i) + deltaqchaud
     155
     156                ! melting + evaporation
     157                zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
     158                ! qbs update
     159                zqbs(i)=max(zqbs(i)-deltaqchaud,0.)
     160             ENDIF
     161
     162             ! now sedimentation scheme with an exact numerical resolution
     163             ! (correct if fall velocity is constant)
     164             zqbsi(i)=zqbs(i)
     165             zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
     166             zqbs(i) = max(zqbs(i),0.)
     167
     168             ! flux update
     169             sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
     170             sedim(i) = max(0.,sedim(i))
    140171     
    141      
    142           ELSE     
     172          ELSE
     173              ! negative celcius temperature     
     174              ! Sublimation scheme   
     175
    143176              zqevti = coef_eva_bs*(1.0-zq(i)/qsi(i))*(sedim(i)**expo_eva_bs) &
    144177                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    … …  
    156189              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime               
    157190              zq(i) = max(0., zq(i))
    158               !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    159               !zqbs(i) = max(0., zqbs(i))
    160191              zt(i) = zt(i)                             &
    161192                + (sedimn(i)-sedim(i))                      &
    162193                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
    163194                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
     195
     196               sedim(i)=sedimn(i)
     197               zqbsi(i)=zqbs(i)
     198
     199               ! now sedimentation scheme with an exact numerical resolution
     200               ! (correct if fall velocity is constant)
     201
     202               zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
     203               zqbs(i) = max(zqbs(i),0.)
     204
     205               ! flux update
     206               sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
     207               sedim(i) = max(0.,sedim(i))
     208
     209
    164210           ENDIF
    165211
    166212
    167           ! sedim update
    168           sedim(i)=sedimn(i)
    169 
    170 
    171         ELSE
    172            sedim(i)=0.
    173         ENDIF ! if sedim > 0
    174 
    175         zqbsi(i)=zqbs(i)
     213           ! if qbs<qbsmin, sublimate or melt and evaporate qbs
     214           ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin
     215           IF (zqbs(i) .LT. qbsmin) THEN                       
     216              zq(i) = zq(i)+zqbs(i)
     217              IF (zt(i) .LT. RTT) THEN
     218                 zt(i) = zt(i) - zqbs(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
     219              ELSE
     220               zt(i) = zt(i) - zqbs(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
     221              ENDIF
     222              zqbs(i)=0.
     223           ENDIF
     224
    176225
    177226    ENDDO ! loop on ngrid
    178227
    179     ! Now sedimention scheme
    180 
    181     ! exact resolution of the conservation equation for qbs with the updated flux from the top (after evap)
    182     ! valid only if the fall velocity is constant
    183 
    184     DO i = 1, ngrid
    185    
    186         zrho(i,k)  = pplay(i,k) / zt(i) / RD
    187         zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
    188         velo(i,k) = fallv_bs
    189         zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
    190         zqbs(i) = max(zqbs(i),0.)
    191         ! flux update
    192         sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
    193         sedim(i) = max(0.,sedim(i))
    194  
    195     ENDDO
    196 
    197 ! old version with bugs
    198 !    DO n = 1, niter_bs
    199 !       DO i = 1, ngrid
    200 !         zrho(i,k)  = pplay(i,k) / zt(i) / RD
    201 !         zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
    202 !         velo(i,k) = fallv_bs
    203 !         dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k)   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
    204 !         precbs   = MIN(MAX(dqsedim,0.0),zqbs(i))
    205 !         zqbs(i) = MAX(zqbs(i)-1.*precbs  , 0.0)
    206 !       ENDDO !loop on ngrid
    207 !    ENDDO ! loop on niter_bs
    208 !    ! add to non sublimated precip
    209 !    DO i=1,ngrid
    210 !       sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    211 !    ENDDO
    212 !
    213228
    214229
    … …  
    218233        bsfl(i,k)=sedim(i)
    219234        dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
    220         dq_bs(i,k) = zq(i) - q(i,k)
     235        dq_bs(i,k) = zq(i) - qv(i,k)
    221236        dtemp_bs(i,k) = zt(i) - temp(i,k)
    222237    ENDDO
    … …  
    235250
    236251end subroutine blowing_snow_sublim_sedim
     252end module lmdz_blowing_snow_sublim_sedim

  • LMDZ6/trunk/libf/phylmd/lmdz_call_blowing_snow.F90

    r4723 r4724  
     1module lmdz_call_blowing_snow
     2
     3contains
     4
    15subroutine call_blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs, &
    26                                  dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
    37
     8use lmdz_blowing_snow_sublim_sedim, only : blowing_snow_sublim_sedim
    49implicit none
    510
    … …  
    3742
    3843end subroutine call_blowing_snow_sublim_sedim
     44
     45end module lmdz_call_blowing_snow

  • LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

    r4722 r4724  
    318318         dser, dt_ds, zsig, zmea
    319319    use phys_output_var_mod, only: tkt, tks, taur, sss
    320     use blowing_snow_ini_mod, only : zeta_bs
     320    use lmdz_blowing_snow_ini, only : zeta_bs
    321321    USE wxios, ONLY: missing_val_xios => missing_val, using_xios
    322322    USE netcdf, only: missing_val_netcdf => nf90_fill_real

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r4722 r4724  
    8080    USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
    8181    USE lmdz_lscp_old, ONLY : fisrtilp
     82    USE lmdz_call_blowing_snow, ONLY : call_blowing_snow_sublim_sedim
    8283    USE lmdz_wake_ini, ONLY : wake_ini
    8384    USE yamada_ini_mod, ONLY : yamada_ini
    … …  
    8586    USE lmdz_thermcell_ini, ONLY : thermcell_ini, iflag_thermals_tenv
    8687    USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
    87     USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs
     88    USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs
    8889    USE lmdz_lscp_ini, ONLY : lscp_ini
    8990    USE lmdz_ratqs_main, ONLY : ratqs_main

  • LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90

    r4673 r4724  
    3535!FC
    3636    USE ioipsl_getin_p_mod, ONLY : getin_p
    37     USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
     37    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
    3838
    3939#ifdef CPP_INLANDSIS

  • LMDZ6/trunk/libf/phylmdiso/lmdz_blowing_snow_ini.F90

    r4723 r4724  
    1 link ../phylmd/blowing_snow_ini_mod.F90
     1link ../phylmd/lmdz_blowing_snow_ini.F90

  • LMDZ6/trunk/libf/phylmdiso/lmdz_blowing_snow_sublim_sedim.F90

    r4723 r4724  
    1 link ../phylmd/blowing_snow_sublim_sedim.F90
     1link ../phylmd/lmdz_blowing_snow_sublim_sedim.F90

  • LMDZ6/trunk/libf/phylmdiso/lmdz_call_blowing_snow.F90

    r4723 r4724  
    1 link ../phylmd/call_blowing_snow_sublim_sedim.F90
     1link ../phylmd/lmdz_call_blowing_snow.F90

  • LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90

    r4723 r4724  
    411411         dser, dt_ds, zsig, zmea
    412412    use phys_output_var_mod, only: tkt, tks, taur, sss
    413     use blowing_snow_ini_mod, only : zeta_bs
     413    use lmdz_blowing_snow_ini, only : zeta_bs
    414414    USE wxios, ONLY: missing_val_xios => missing_val, using_xios
    415415    use netcdf, only: nf90_fill_real

  • LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

    r4723 r4724  
    8080    USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
    8181    USE lmdz_lscp_old, ONLY : fisrtilp
     82    USE lmdz_call_blowing_snow, ONLY : call_blowing_snow_sublim_sedim
    8283    USE lmdz_wake_ini, ONLY : wake_ini
    8384    USE yamada_ini_mod, ONLY : yamada_ini
    … …  
    8586    USE lmdz_thermcell_ini, ONLY : thermcell_ini, iflag_thermals_tenv
    8687    USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
    87     USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs
     88    USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs
    8889    USE lmdz_lscp_ini, ONLY : lscp_ini
    8990    USE lmdz_ratqs_main, ONLY : ratqs_main

  • LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90

    r4673 r4724  
    5050!FC
    5151    USE ioipsl_getin_p_mod, ONLY : getin_p
    52     USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
     52    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
    5353
    5454#ifdef CPP_INLANDSIS
Note: See TracChangeset for help on using the changeset viewer.