Changeset 5489 for LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90

Timestamp:

Jan 17, 2025, 6:16:25 PM (13 days ago)

Author:

aborella

Message:

Merge with trunk

Location:

LMDZ6/branches/contrails

Files:

• . (modified) (1 prop)
• libf/phylmd/physiq_mod.F90 (modified) (20 diffs)

LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90

@@ -39 +39 @@
     USE ioipsl_getin_p_mod, ONLY : getin_p
     USE indice_sol_mod
-    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, addPhase
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, addPhase, ivap, iliq, isol, ibs, icf, irvc, ircont
     USE strings_mod,  ONLY: strIdx
     USE iophy
@@ -78 +78 @@
     USE lmdz_lscp, ONLY : lscp
     USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
-    USE lmdz_lscp_old, ONLY : fisrtilp
+    USE lmdz_lscp_old, ONLY : fisrtilp, fisrtilp_first
     USE lmdz_call_blowing_snow, ONLY : call_blowing_snow_sublim_sedim
     USE lmdz_wake_ini, ONLY : wake_ini
@@ -248 +248 @@
        cldh, cldl,cldm, cldq, cldt,      &
        JrNt,                             &
-       dthmin, evap, snowerosion,fder, plcl, plfc,   &
+       dthmin, evap, snowerosion, icesub_lic, fder, plcl, plfc,   &
        prw, prlw, prsw, prbsw, water_budget,         &
        s_lcl, s_pblh, s_pblt, s_therm,   &
@@ -376 +376 @@
        USE phys_output_write_spl_mod, ONLY: phys_output_write_spl
        USE phytracr_spl_mod, ONLY: phytracr_spl_out_init, phytracr_spl
+       USE s2s, ONLY : s2s_initialize
     IMPLICIT NONE
     !>======================================================================
@@ -512 +513 @@
     !======================================================================
     !
-    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
-    INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, irvc, ircont
-!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, irvc, ircont)
-    !
     !
     ! Variables argument:
@@ -1021 +1018 @@
 
     REAL picefra(klon,klev)
-    REAL zrel_oro(klon)
+    REAL nm_oro(klon)
     !IM cf. AM 081204 END
     !
@@ -1096 +1093 @@
     CHARACTER*80 abort_message
     LOGICAL, SAVE ::  ok_sync, ok_sync_omp
-    !$OMP THREADPRIVATE(ok_sync)
+    !$OMP THREADPRIVATE(ok_sync,ok_sync_omp)
     REAL date0
 
@@ -1106 +1103 @@
     REAL ztsol(klon)
     REAL q2m(klon,nbsrf)  ! humidite a 2m
-    REAL fsnowerosion(klon,nbsrf) ! blowing snow flux at surface
     REAL qbsfra  ! blowing snow fraction
     !IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
@@ -1270 +1266 @@
     ! Subgrid scale wind :
     ! Need to be allocatable/save because the number of bin is not known (provided by surf_wind_ini)
-    integer, save :: nsrfwnd=1
+    integer, save :: nsurfwind=1
     real, dimension(:,:), allocatable, save :: surf_wind_value, surf_wind_proba ! module and probability of sugrdi wind wind sample
-    !$OMP THREADPRIVATE(nsrfwnd,surf_wind_value, surf_wind_proba)
+    !$OMP THREADPRIVATE(nsurfwind,surf_wind_value, surf_wind_proba)
     
 
@@ -1352 +1348 @@
 
     IF (first) THEN
-       ivap = strIdx(tracers(:)%name, addPhase('H2O', 'g'))
-       iliq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
-       isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
-       ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
-       icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
-       irvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
-       ircont = strIdx(tracers(:)%name, addPhase('H2O', 'a'))
+        
+        CALL s2s_initialize     ! initialization of source to source tools
+        
 !       CALL init_etat0_limit_unstruct
 !       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
@@ -1841 +1833 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        CALL surf_wind_ini(klon,lunout)
-       CALL getin_p('nsrfwnd',nsrfwnd)
-       allocate(surf_wind_value(klon,nsrfwnd),surf_wind_proba(klon,nsrfwnd))
+       CALL getin_p('nsurfwind',nsurfwind)
+       allocate(surf_wind_value(klon,nsurfwind),surf_wind_proba(klon,nsurfwind))
     
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   CALL iophys_ini(pdtphys,nsurfwind) ! replay automatic include  ! replay automatic include
        CALL wake_ini(rg,rd,rv,prt_level)
        CALL yamada_ini(klon,lunout,prt_level)
@@ -2917 +2910 @@
             cdragh,    cdragm,  u1,    v1,            &
             beta_aridity, &
-                                !albedo SB >>>
-                                ! albsol1,   albsol2,   sens,    evap,      &
-            albsol_dir,   albsol_dif,   sens,    evap, snowerosion, &
-                                !albedo SB <<<
+            albsol_dir,   albsol_dif,   sens,    evap, snowerosion, icesub_lic, &
             albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
             zxtsol,    zxfluxlat, zt2m,    qsat2m,  zn2mout, &
@@ -3736 +3726 @@
              !  poches, la tendance moyenne associ\'ee doit etre
              !  multipliee par la fraction surfacique qu'ils couvrent.
+             IF (mod(iflag_pbl_split/10,10) == 1) THEN
+                ! On tient compte du splitting pour modifier les profils deltatq/T des poches
+                DO k=1,klev
+                   DO i=1,klon
+                      d_deltat_the(i,k) = - d_t_ajs(i,k)
+                      d_deltaq_the(i,k) = - d_q_ajs(i,k)
+                   ENDDO
+                ENDDO
+             ELSE
+                d_deltat_the(:,:) = 0.
+                d_deltaq_the(:,:) = 0.
+             ENDIF
+
              DO k=1,klev
                 DO i=1,klon
-                   !
-                   d_deltat_the(i,k) = - d_t_ajs(i,k)
-                   d_deltaq_the(i,k) = - d_q_ajs(i,k)
-                   !
                    d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
                    d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
                    d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
                    d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
-                   !
                 ENDDO
              ENDDO
@@ -3837 +3835 @@
     !===================================================================
     ! Computation of subrgid scale near-surface wind distribution
-    call surf_wind(klon,nsrfwnd,u10m,v10m,wake_s,wake_Cstar,ustar,wstar,surf_wind_value,surf_wind_proba)
+    call surf_wind(klon,nsurfwind,u10m,v10m,wake_s,wake_Cstar,ustar,wstar,surf_wind_value,surf_wind_proba)
 
     !===================================================================
@@ -3924 +3922 @@
 
     ELSE
-
+    
+    CALL fisrtilp_first(klon, klev, phys_tstep, pfrac_impa, pfrac_nucl, pfrac_1nucl)
     CALL fisrtilp(klon,klev,phys_tstep,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs,sigma_qtherm, &
@@ -4859 +4858 @@
     ! a l'echelle sous-maille:
     !
+
+    ! calculation of nm_oro
+    DO i=1,klon
+          ! nm_oro is a proxy for the number of subgrid scale mountains
+          ! -> condition on nm_oro can deactivate the lifting on tilted planar terrains
+          !    such as ice sheets (work by V. Wiener)
+          ! in such a case, the SSO scheme should activate only where nm_oro>0 i.e. by setting
+          ! nm_oro_t=0.
+          nm_oro(i)=zsig(i)*sqrt(cell_area(i)*(pctsrf(i,is_ter)+pctsrf(i,is_lic)))/(4.*MAX(zstd(i),1.e-8))-1.
+    ENDDO
+
     IF (prt_level .GE.10) THEN
        print *,' call orography ? ', ok_orodr
@@ -4869 +4879 @@
        DO i=1,klon
           itest(i)=0
-          zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
-          !zrel_oro: relative mountain height wrt relief explained by mean slope
-          ! -> condition on zrel_oro can deactivate the drag on tilted planar terrains
-          !    such as ice sheets (work by V. Wiener)
           ! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to
           ! earn computation time but they are not physical.
-          IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
+          IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(nm_oro(i).GT.nm_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -4924 +4930 @@
        DO i=1,klon
           itest(i)=0
-          !zrel_oro: relative mountain height wrt relief explained by mean slope
-          ! -> condition on zrel_oro can deactivate the lifting on tilted planar terrains
-          !    such as ice sheets (work by V. Wiener)
-          zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
-          IF (((zpic(i)-zmea(i)).GT.zpmm_orolf_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
+          IF (((zpic(i)-zmea(i)).GT.zpmm_orolf_t).AND.(nm_oro(i).GT.nm_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -5169 +5171 @@
 ! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
 ! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
-          IF ((zstd(i).GT.1.0) .AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
+          IF ((zstd(i).GT.1.0) .AND.(nm_oro(i).GT.nm_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -5181 +5183 @@
        DO i=1,klon
           itest(i)=0
-        IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
+        IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(nm_oro(i).GT.nm_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1