Changeset 4724 for LMDZ6/trunk/libf

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:

• phylmd/lmdz_blowing_snow_ini.F90 (moved) (moved from LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90) (3 diffs)
• phylmd/lmdz_blowing_snow_sublim_sedim.F90 (moved) (moved from LMDZ6/trunk/libf/phylmd/blowing_snow_sublim_sedim.F90) (9 diffs)
• phylmd/lmdz_call_blowing_snow.F90 (moved) (moved from LMDZ6/trunk/libf/phylmd/call_blowing_snow_sublim_sedim.F90) (2 diffs)
• phylmd/pbl_surface_mod.F90 (modified) (1 diff)
• phylmd/physiq_mod.F90 (modified) (2 diffs)
• phylmd/surf_landice_mod.F90 (modified) (1 diff)
• phylmdiso/lmdz_blowing_snow_ini.F90 (moved) (moved from LMDZ6/trunk/libf/phylmdiso/blowing_snow_ini_mod.F90) (1 diff)
• phylmdiso/lmdz_blowing_snow_sublim_sedim.F90 (moved) (moved from LMDZ6/trunk/libf/phylmdiso/blowing_snow_sublim_sedim.F90) (1 diff)
• phylmdiso/lmdz_call_blowing_snow.F90 (moved) (moved from LMDZ6/trunk/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90) (1 diff)
• phylmdiso/pbl_surface_mod.F90 (modified) (1 diff)
• phylmdiso/physiq_mod.F90 (modified) (2 diffs)
• phylmdiso/surf_landice_mod.F90 (modified) (1 diff)

LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90

@@ -1 @@
-module blowing_snow_ini_mod
+module lmdz_blowing_snow_ini
 
 implicit none
@@ -15 +15 @@
    !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs)
    !$OMP THREADPRIVATE(iflag_saltation_bs)
+
+    real tbsmelt                  ! parameter to calculate melting fraction of BS sedimentation
+    parameter (tbsmelt=278.15)
+
+    real taumeltbs0               ! Melting time scale of blowing snow at 273.15K
+    parameter (taumeltbs0=1800.0)
+
+    real qbsmin                   ! Minimum blowing snow specific content
+    parameter (qbsmin=1.E-10)
+
 
       contains
@@ -70 +80 @@
       end subroutine blowing_snow_ini
 
-end module blowing_snow_ini_mod
+end module lmdz_blowing_snow_ini

LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90

@@ -1 @@
-subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
+module lmdz_blowing_snow_sublim_sedim
+
+contains
+subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
 
 !==============================================================================
@@ -8 +11 @@
 
 
-use blowing_snow_ini_mod, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs
-use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2
+use lmdz_blowing_snow_ini, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs, qbsmin
+use lmdz_blowing_snow_ini, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, tbsmelt, taumeltbs0
 USE lmdz_lscp_tools, only : calc_qsat_ecmwf
 
@@ -25 +28 @@
 real, intent(in)                        :: dtime
 real, intent(in), dimension(ngrid,nlay) :: temp
-real, intent(in), dimension(ngrid,nlay) :: q
+real, intent(in), dimension(ngrid,nlay) :: qv
 real, intent(in), dimension(ngrid,nlay) :: qbs
 real, intent(in), dimension(ngrid,nlay) :: pplay
@@ -49 +52 @@
 integer                                  :: k,i,n
 real                                     :: zqev0, zqevi, zqevti, zcpair, zcpeau, dqbsmelt
-real                                     :: dqsedim,precbs
+real                                     :: dqsedim,precbs, deltaqchaud, zmelt, taumeltbs
+real                                     :: maxdeltaqchaud
 real, dimension(ngrid)                   :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn
 real, dimension(ngrid)                   :: zqbsi,zmqc, zmair, zdz
@@ -73 +77 @@
     DO i=1,ngrid
         zt(i)=temp(i,k)
-        zq(i)=q(i,k)
+        zq(i)=qv(i,k)
         zqbs(i)=qbs(i,k)
     ENDDO
@@ -111 +115 @@
     
     DO i = 1, ngrid
-       ! if sedimentation:
-       IF (sedim(i) .GT. 0.) THEN
+
+       zrho(i,k)  = pplay(i,k) / zt(i) / RD
+       zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
+       ! BS fall velocity
+       velo(i,k) = fallv_bs
+
 
           IF (zt(i) .GT. RTT) THEN
+             
              ! if positive celcius temperature, we assume
-             ! that all the blowing snow melt and evaporate
-              zqevti=sedim(i)*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG     
-             ! we ensure that the whole mesh does not exceed saturation wrt liquid
-              zqev0 = MAX(0.0, qsl(i)-zq(i))
-              zqevi = MIN(zqev0,zqevti)
-              ! New solid precipitation fluxes
-              sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
-              ! vapor, temperature, precip fluxes update
-              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-              zq(i) = max(0., zq(i))
-              !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-              !zqbs(i) = max(0., zqbs(i))
-              ! melting
-              zt(i) = zt(i)                             &
-                + (sedimn(i)-sedim(i))                      &
-                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
-                * RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
-              ! evaporation
-              zt(i) = zt(i)                             &
-                + (sedimn(i)-sedim(i))                      &
-                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
-                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))      
+             ! that part of the the blowing snow flux melts and evaporates
+             
+             ! vapor, bs, temperature, precip fluxes update
+             zmelt = ((zt(i)-RTT)/(tbsmelt-RTT))
+             zmelt = MIN(MAX(zmelt,0.),1.)
+             sedimn(i)=sedim(i)*zmelt
+             deltaqchaud=-(sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+             ! max evap such as celcius temperature cannot become negative
+             maxdeltaqchaud= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
+
+             deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud)
+             zq(i) = zq(i) + deltaqchaud  
+              
+             ! melting + evaporation
+             zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+
+             sedim(i)=sedimn(i)
+
+             ! if temperature still positive, we assume that part of the blowing snow 
+             ! already present in the mesh melts and evaporates with a typical time 
+             ! constant between taumeltbs0 and 0
+ 
+             IF ( zt(i) .GT. RTT ) THEN
+                taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT)))
+                deltaqchaud=min(zqbs(i),zqbs(i)/taumeltbs*dtime)
+                maxdeltaqchaud= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
+                deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud)
+                zq(i) = zq(i) + deltaqchaud
+
+                ! melting + evaporation
+                zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+                ! qbs update
+                zqbs(i)=max(zqbs(i)-deltaqchaud,0.)
+             ENDIF
+
+             ! now sedimentation scheme with an exact numerical resolution 
+             ! (correct if fall velocity is constant)
+             zqbsi(i)=zqbs(i)
+             zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
+             zqbs(i) = max(zqbs(i),0.)
+
+             ! flux update
+             sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
+             sedim(i) = max(0.,sedim(i))
       
-      
-          ELSE     
+          ELSE 
+              ! negative celcius temperature     
+              ! Sublimation scheme   
+
               zqevti = coef_eva_bs*(1.0-zq(i)/qsi(i))*(sedim(i)**expo_eva_bs) &
                     *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
@@ -156 +189 @@
               zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime                
               zq(i) = max(0., zq(i))
-              !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-              !zqbs(i) = max(0., zqbs(i))
               zt(i) = zt(i)                             &
                 + (sedimn(i)-sedim(i))                      &
                 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
                 * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+
+               sedim(i)=sedimn(i)
+               zqbsi(i)=zqbs(i)
+
+               ! now sedimentation scheme with an exact numerical resolution 
+               ! (correct if fall velocity is constant)
+
+               zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
+               zqbs(i) = max(zqbs(i),0.)
+
+               ! flux update
+               sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
+               sedim(i) = max(0.,sedim(i))
+
+
            ENDIF
 
 
-          ! sedim update
-          sedim(i)=sedimn(i)
-
-
-        ELSE
-           sedim(i)=0.
-        ENDIF ! if sedim > 0
-
-        zqbsi(i)=zqbs(i)
+           ! if qbs<qbsmin, sublimate or melt and evaporate qbs
+           ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin
+           IF (zqbs(i) .LT. qbsmin) THEN                       
+              zq(i) = zq(i)+zqbs(i)
+              IF (zt(i) .LT. RTT) THEN
+                 zt(i) = zt(i) - zqbs(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+              ELSE
+               zt(i) = zt(i) - zqbs(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+              ENDIF
+              zqbs(i)=0.
+           ENDIF
+
 
     ENDDO ! loop on ngrid
 
-    ! Now sedimention scheme 
-
-    ! exact resolution of the conservation equation for qbs with the updated flux from the top (after evap)
-    ! valid only if the fall velocity is constant
-
-    DO i = 1, ngrid
-   
-        zrho(i,k)  = pplay(i,k) / zt(i) / RD
-        zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
-        velo(i,k) = fallv_bs
-        zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
-        zqbs(i) = max(zqbs(i),0.)
-        ! flux update
-        sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
-        sedim(i) = max(0.,sedim(i))
- 
-    ENDDO
-
-! old version with bugs
-!    DO n = 1, niter_bs
-!       DO i = 1, ngrid
-!         zrho(i,k)  = pplay(i,k) / zt(i) / RD
-!         zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
-!         velo(i,k) = fallv_bs 
-!         dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k)   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
-!         precbs   = MIN(MAX(dqsedim,0.0),zqbs(i))
-!         zqbs(i) = MAX(zqbs(i)-1.*precbs  , 0.0)
-!       ENDDO !loop on ngrid
-!    ENDDO ! loop on niter_bs
-!    ! add to non sublimated precip
-!    DO i=1,ngrid
-!       sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
-!    ENDDO
-!
 
 
@@ -218 +233 @@
         bsfl(i,k)=sedim(i)
         dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
-        dq_bs(i,k) = zq(i) - q(i,k)
+        dq_bs(i,k) = zq(i) - qv(i,k)
         dtemp_bs(i,k) = zt(i) - temp(i,k)
     ENDDO
@@ -235 +250 @@
 
 end subroutine blowing_snow_sublim_sedim
+end module lmdz_blowing_snow_sublim_sedim

LMDZ6/trunk/libf/phylmd/lmdz_call_blowing_snow.F90

@@ +1 @@
+module lmdz_call_blowing_snow
+
+contains
+
 subroutine call_blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs, &
                                   dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
 
+use lmdz_blowing_snow_sublim_sedim, only : blowing_snow_sublim_sedim
 implicit none
 
@@ -37 +42 @@
 
 end subroutine call_blowing_snow_sublim_sedim
+
+end module lmdz_call_blowing_snow

LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

@@ -318 +318 @@
          dser, dt_ds, zsig, zmea
     use phys_output_var_mod, only: tkt, tks, taur, sss
-    use blowing_snow_ini_mod, only : zeta_bs
+    use lmdz_blowing_snow_ini, only : zeta_bs
     USE wxios, ONLY: missing_val_xios => missing_val, using_xios
     USE netcdf, only: missing_val_netcdf => nf90_fill_real

LMDZ6/trunk/libf/phylmd/physiq_mod.F90

@@ -80 +80 @@
     USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
     USE lmdz_lscp_old, ONLY : fisrtilp
+    USE lmdz_call_blowing_snow, ONLY : call_blowing_snow_sublim_sedim
     USE lmdz_wake_ini, ONLY : wake_ini
     USE yamada_ini_mod, ONLY : yamada_ini
@@ -85 +86 @@
     USE lmdz_thermcell_ini, ONLY : thermcell_ini, iflag_thermals_tenv
     USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
-    USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs 
+    USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs 
     USE lmdz_lscp_ini, ONLY : lscp_ini
     USE lmdz_ratqs_main, ONLY : ratqs_main

LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90

@@ -35 +35 @@
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
+    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
 
 #ifdef CPP_INLANDSIS

LMDZ6/trunk/libf/phylmdiso/lmdz_blowing_snow_ini.F90

@@ -1 @@
-link ../phylmd/blowing_snow_ini_mod.F90
+link ../phylmd/lmdz_blowing_snow_ini.F90

LMDZ6/trunk/libf/phylmdiso/lmdz_blowing_snow_sublim_sedim.F90

@@ -1 @@
-link ../phylmd/blowing_snow_sublim_sedim.F90
+link ../phylmd/lmdz_blowing_snow_sublim_sedim.F90

LMDZ6/trunk/libf/phylmdiso/lmdz_call_blowing_snow.F90

@@ -1 @@
-link ../phylmd/call_blowing_snow_sublim_sedim.F90
+link ../phylmd/lmdz_call_blowing_snow.F90

LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90

@@ -411 +411 @@
          dser, dt_ds, zsig, zmea
     use phys_output_var_mod, only: tkt, tks, taur, sss
-    use blowing_snow_ini_mod, only : zeta_bs
+    use lmdz_blowing_snow_ini, only : zeta_bs
     USE wxios, ONLY: missing_val_xios => missing_val, using_xios
     use netcdf, only: nf90_fill_real

LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

@@ -80 +80 @@
     USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
     USE lmdz_lscp_old, ONLY : fisrtilp
+    USE lmdz_call_blowing_snow, ONLY : call_blowing_snow_sublim_sedim
     USE lmdz_wake_ini, ONLY : wake_ini
     USE yamada_ini_mod, ONLY : yamada_ini
@@ -85 +86 @@
     USE lmdz_thermcell_ini, ONLY : thermcell_ini, iflag_thermals_tenv
     USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
-    USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs 
+    USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs 
     USE lmdz_lscp_ini, ONLY : lscp_ini
     USE lmdz_ratqs_main, ONLY : ratqs_main

LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90

@@ -50 +50 @@
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
+    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
 
 #ifdef CPP_INLANDSIS