Ignore:
Timestamp:
May 23, 2024, 12:02:33 PM (8 months ago)
Author:
aborella
Message:

New version of condensation and ice supersaturation in LSCP.
Multiple changes troughout the code (in particular, two new water phase tracers).

Location:
LMDZ6/branches/cirrus/libf/phylmd
Files:
1 added
1 deleted
13 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/cirrus/libf/phylmd/clesphys.h

    r4773 r4951  
    9494       LOGICAL :: ok_airs
    9595       INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo, NSW, iflag_albedo
    96        LOGICAL :: ok_ice_sursat, ok_plane_h2o, ok_plane_contrail
     96       LOGICAL :: ok_ice_supersat, ok_plane_h2o, ok_plane_contrail
    9797       LOGICAL :: ok_chlorophyll
    9898       LOGICAL :: ok_strato
    … …  
    154154     &     , ok_lic_melt, ok_lic_cond, aer_type                         &
    155155     &     , iflag_rrtm, ok_strato,ok_hines, ok_qch4                    &
    156      &     , iflag_ice_thermo, ok_ice_sursat                            &
     156     &     , iflag_ice_thermo, ok_ice_supersat                          &
    157157     &     , ok_plane_h2o, ok_plane_contrail                            &
    158158     &     , ok_gwd_rando, NSW, iflag_albedo                            &

  • LMDZ6/branches/cirrus/libf/phylmd/conf_phys_m.F90

    r4843 r4951  
    173173    INTEGER,SAVE :: iflag_clw_omp
    174174    INTEGER,SAVE :: iflag_ice_thermo_omp
    175     LOGICAL,SAVE :: ok_ice_sursat_omp
     175    LOGICAL,SAVE :: ok_ice_supersat_omp
    176176    LOGICAL,SAVE :: ok_plane_h2o_omp, ok_plane_contrail_omp
    177177    REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
    … …  
    13471347
    13481348    !
    1349     !Config Key  = ok_ice_sursat
    1350     !Config Desc =
    1351     !Config Def  = 0
    1352     !Config Help =
    1353     !
    1354     ok_ice_sursat_omp = 0
    1355     CALL getin('ok_ice_sursat',ok_ice_sursat_omp)
     1349    !Config Key  = ok_ice_supersat
     1350    !Config Desc = include ice supersaturation for cold clouds
     1351    !Config Def  = .FALSE.
     1352    !Config Help =
     1353    !
     1354    ok_ice_supersat_omp = .FALSE.
     1355    CALL getin('ok_ice_supersat',ok_ice_supersat_omp)
    13561356
    13571357    !Config Key  = ok_plane_h2o
    1358     !Config Desc =
    1359     !Config Def  = 0
     1358    !Config Desc = include H2O emissions from aviation
     1359    !Config Def  = .FALSE.
    13601360    !Config Help =
    13611361    !
    … …  
    13641364
    13651365    !Config Key  = ok_plane_contrail
    1366     !Config Desc =
    1367     !Config Def  = 0
     1366    !Config Desc = include the formation of contrail cirrus clouds
     1367    !Config Def  = .FALSE.
    13681368    !Config Help =
    13691369    !
    … …  
    23452345    rad_chau2 = rad_chau2_omp
    23462346    iflag_ice_thermo = iflag_ice_thermo_omp
    2347     ok_ice_sursat = ok_ice_sursat_omp
     2347    ok_ice_supersat = ok_ice_supersat_omp
    23482348    ok_plane_h2o = ok_plane_h2o_omp
    23492349    ok_plane_contrail = ok_plane_contrail_omp
    … …  
    27702770    WRITE(lunout,*) ' rad_chau2 = ',rad_chau2
    27712771    WRITE(lunout,*) ' iflag_ice_thermo = ',iflag_ice_thermo
    2772     WRITE(lunout,*) ' ok_ice_sursat = ',ok_ice_sursat
     2772    WRITE(lunout,*) ' ok_ice_supersat = ',ok_ice_supersat
    27732773    WRITE(lunout,*) ' ok_plane_h2o = ',ok_plane_h2o
    27742774    WRITE(lunout,*) ' ok_plane_contrail = ',ok_plane_contrail

  • LMDZ6/branches/cirrus/libf/phylmd/create_etat0_unstruct_mod.F90

    r4856 r4951  
    258258    prw_ancien = 0.
    259259    agesno = 0.
     260   
     261    ! LSCP condensation and ice supersaturation
     262    cf_ancien = 0.
     263    rvc_ancien = 0.
    260264
    261265    wake_delta_pbl_TKE(:,:,:)=0
    … …  
    310314    END IF
    311315    ratqs_inter_ = 0.002
    312     rneb_ancien = 0.
    313316
    314317    CALL gather_omp(cell_area,cell_area_mpi)

  • LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90

    r4915 r4951  
    88SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
    99     paprs,pplay,temp,qt,ptconv,ratqs,                  &
    10      d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri,  &
     10     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
    1111     pfraclr,pfracld,                                   &
    1212     radocond, radicefrac, rain, snow,                  &
    … …  
    1414     prfl, psfl, rhcl, qta, fraca,                     &
    1515     tv, pspsk, tla, thl, iflag_cld_th,             &
    16      iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
    17      distcltop,temp_cltop,                              &
    18      qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
    19      Tcontr, qcontr, qcontr2, fcontrN, fcontrP,         &
     16     iflag_ice_thermo, distcltop, temp_cltop, cell_area,&
     17     cf_seri, rvc_seri, u_seri, v_seri, pbl_eps,        &
     18     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
     19     ratio_qi_qtot, dcf_sub, dcf_con, dcf_mix,          &
     20     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
     21     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
     22     Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi,&
     23     dqi_avi, dqvc_avi, flight_dist, flight_h2o,        &
    2024     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
    2125     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
    … …  
    98102USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat
    99103USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
    100 USE ice_sursat_mod, ONLY : ice_sursat
     104USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
    101105USE lmdz_lscp_poprecip, ONLY : poprecip_precld, poprecip_postcld
    102106
    103107! Use du module lmdz_lscp_ini contenant les constantes
    104 USE lmdz_lscp_ini, ONLY : prt_level, lunout
     108USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
    105109USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
    106110USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
    … …  
    111115USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
    112116USE lmdz_lscp_ini, ONLY : ok_poprecip
     117USE lmdz_lscp_ini, ONLY : ok_external_lognormal, ok_ice_supersat, ok_unadjusted_clouds
    113118
    114119IMPLICIT NONE
    … …  
    132137  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
    133138                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
    134   LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated
    135 
    136139  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
    137140
    … …  
    150153  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs        ! function of pressure that sets the large-scale
    151154
    152 
    153   ! Input sursaturation en glace
    154   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rneb_seri        ! fraction nuageuse en memoire
     155  ! INPUT/OUTPUT condensation and ice supersaturation
     156  !--------------------------------------------------
     157  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
     158  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratio_qi_qtot    ! solid specific water to total specific water ratio [-]
     159  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
     160  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
     161  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
     162  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: pbl_eps          ! TKE dissipation [?]
     163  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
     164
     165  ! INPUT/OUTPUT aviation
     166  !--------------------------------------------------
     167  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
     168  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
    155169 
    156170  ! OUTPUT variables
    … …  
    170184  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
    171185  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
    172   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl            ! saturation specific humidity wrt liquid [kg/kg]
    173   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsats            ! saturation specific humidity wrt ice [kg/kg] 
    174186  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
    175187  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
    … …  
    183195  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl        ! scavenging fraction due tu nucleation [-]
    184196 
    185   ! for supersaturation and contrails parameterisation
    186  
    187   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qclr             ! specific total water content in clear sky region [kg/kg]
    188   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld             ! specific total water content in cloudy region [kg/kg]
    189   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qss              ! specific total water content in supersat region [kg/kg]
    190   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qvc              ! specific vapor content in clouds [kg/kg]
    191   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebclr          ! mesh fraction of clear sky [-]   
    192   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebss           ! mesh fraction of ISSR [-] 
    193   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_ss         ! coefficient governing the ice nucleation RHi threshold [-]     
    194   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr           ! threshold temperature for contrail formation [K]
    195   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr           ! threshold humidity for contrail formation [kg/kg]
    196   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2          ! // (2nd expression more consistent with LMDZ expression of q)         
    197   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN          ! fraction of grid favourable to non-persistent contrails
    198   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP          ! fraction of grid favourable to persistent contrails
    199   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      ! mean saturation deficit in thermals
    200   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     ! mean saturation deficit in environment
    201   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  ! std of saturation deficit in thermals
    202   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv ! std of saturation deficit in environment
     197  ! for condensation and ice supersaturation
     198
     199  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
     200  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
     201  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
     202  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
     203  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
     204  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
     205  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
     206  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
     207  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
     208  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
     209  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
     210  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
     211  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
     212  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
     213  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
     214  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
     215  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
     216  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
     217  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
     218
     219  ! for contrails and aviation
     220
     221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr         !--threshold temperature for contrail formation [K]
     222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr         !--threshold humidity for contrail formation [kg/kg]
     223  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2        !--// (2nd expression more consistent with LMDZ expression of q)
     224  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN        !--fraction of grid favourable to non-persistent contrails
     225  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP        !--fraction of grid favourable to persistent contrails
     226  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
     227  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
     228  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]
     229
    203230
    204231  ! for POPRECIP
    … …  
    217244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
    218245  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
     246
     247  ! for thermals
     248
     249  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
     250  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
     251  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
     252  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
    219253
    220254
    … …  
    282316  REAL :: effective_zneb
    283317  REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D
     318 
     319  ! for condensation and ice supersaturation
     320  REAL, DIMENSION(klon) :: qvc, shear
     321  REAL :: delta_z
     322  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
     323  ! Constants used for calculating ratios that are advected (using a parent-child
     324  ! formalism). This is not done in the dynamical core because at this moment,
     325  ! only isotopes can use this parent-child formalism. Note that the two constants
     326  ! are the same as the one use in the dynamical core, being also defined in
     327  ! dyn3d_common/infotrac.F90
     328  REAL :: min_qParent, min_ratio
    284329
    285330
    … …  
    363408temp_cltop1D(:) = 0.0
    364409ztupnew(:)=0.0
    365 !--ice supersaturation
    366 gamma_ss(:,:) = 1.0
    367 qss(:,:) = 0.0
    368 rnebss(:,:) = 0.0
    369 Tcontr(:,:) = missing_val
    370 qcontr(:,:) = missing_val
    371 qcontr2(:,:) = missing_val
    372 fcontrN(:,:) = 0.0
    373 fcontrP(:,:) = 0.0
     410
    374411distcltop(:,:)=0.
    375412temp_cltop(:,:)=0.
     413
     414!--Ice supersaturation
     415gamma_cond(:,:) = 1.
     416qissr(:,:)      = 0.
     417issrfra(:,:)    = 0.
     418dcf_sub(:,:)    = 0.
     419dcf_con(:,:)    = 0.
     420dcf_mix(:,:)    = 0.
     421dqi_adj(:,:)    = 0.
     422dqi_sub(:,:)    = 0.
     423dqi_con(:,:)    = 0.
     424dqi_mix(:,:)    = 0.
     425dqvc_adj(:,:)   = 0.
     426dqvc_sub(:,:)   = 0.
     427dqvc_con(:,:)   = 0.
     428dqvc_mix(:,:)   = 0.
     429fcontrN(:,:)    = 0.
     430fcontrP(:,:)    = 0.
     431Tcontr(:,:)     = missing_val
     432qcontr(:,:)     = missing_val
     433qcontr2(:,:)    = missing_val
     434dcf_avi(:,:)    = 0.
     435dqi_avi(:,:)    = 0.
     436dqvc_avi(:,:)   = 0.
     437qvc(:)          = 0.
     438shear(:)        = 0.
     439min_qParent     = 1.e-30
     440min_ratio       = 1.e-16
     441
    376442!-- poprecip
    377443qraindiag(:,:)= 0.
    … …  
    759825                    rneblsvol(i,k)=ctot_vol(i,k)
    760826                    zqn(i)=qcloud(i)
     827                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
     828                    qvc(i) = rneb(i,k) * zqs(i)
    761829                ENDDO
    762830
    … …  
    794862        DO iter=1,iflag_fisrtilp_qsat+1
    795863
     864                keepgoing(:) = .FALSE.
     865
    796866                DO i=1,klon
    797867
    798868                    ! keepgoing = .true. while convergence is not satisfied
    799                     keepgoing(i)=ABS(DT(i)).GT.DDT0
    800 
    801                     IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
     869
     870                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
    802871
    803872                        ! if not convergence:
     873                        ! we calculate a new iteration
     874                        keepgoing(i) = .TRUE.
    804875
    805876                        ! P2.2.1> cloud fraction and condensed water mass calculation
    … …  
    833904                  CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:))
    834905
     906
     907                  !--AB Activates a condensation scheme that allows for
     908                  !--ice supersaturation and contrails evolution from aviation
     909                  IF (ok_ice_supersat) THEN
     910
     911                    !--Calculate the shear value (input for condensation and ice supersat)
     912                    DO i = 1, klon
     913                      !--Cell thickness [m]
     914                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
     915                      IF ( iftop ) THEN
     916                        ! top
     917                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
     918                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
     919                      ELSEIF ( k .EQ. 1 ) THEN
     920                        ! surface
     921                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
     922                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
     923                      ELSE
     924                        ! other layers
     925                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
     926                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
     927                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
     928                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
     929                      ENDIF
     930                    ENDDO
     931
     932                    !---------------------------------------------
     933                    !--   CONDENSATION AND ICE SUPERSATURATION  --
     934                    !---------------------------------------------
     935
     936                    CALL condensation_ice_supersat( &
     937                        klon, dtime, missing_val, &
     938                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
     939                        cf_seri(:,k), rvc_seri(:,k), ratio_qi_qtot(:,k), &
     940                        shear(:), pbl_eps(:,k), cell_area(:), &
     941                        Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), keepgoing(:), &
     942                        rneb(:,k), zqn(:), qvc(:), issrfra(:,k), qissr(:,k), &
     943                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
     944                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
     945                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
     946                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
     947                        flight_dist(:,k), flight_h2o(:,k), &
     948                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
     949
     950
     951                  !--If .TRUE., calls an externalised version of the generalised
     952                  !--lognormal condensation scheme (Bony and Emanuel 2001)
     953                  !--Numerically, convergence is conserved with this option
     954                  !--The objective is to simplify LSCP
     955                  ELSEIF ( ok_external_lognormal ) THEN
     956                         
     957                   CALL condensation_lognormal( &
     958                       klon, Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), &
     959                       keepgoing(:), rneb(:,k), zqn(:), qvc(:))
     960
     961
     962                 ELSE !--Generalised lognormal (Bony and Emanuel 2001)
     963
    835964                  DO i=1,klon !todoan : check if loop in i is needed
    836965
    837                       IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
     966                      IF (keepgoing(i)) THEN
    838967
    839968                        zpdf_sig(i)=ratqs(i,k)*zq(i)
    … …  
    849978                        zpdf_e2(i)=1.-erf(zpdf_e2(i))
    850979
    851                         IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN
    852 
    853980                          IF (zpdf_e1(i).LT.1.e-10) THEN
    854981                              rneb(i,k)=0.
    855                               zqn(i)=gammasat(i)*zqs(i)
     982                              zqn(i)=zqs(i)
     983                              !--AB grid-mean vapor in the cloud - we assume saturation adjustment
     984                              qvc(i) = 0.
    856985                          ELSE
    857986                              rneb(i,k)=0.5*zpdf_e1(i)
    858987                              zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
     988                              !--AB grid-mean vapor in the cloud - we assume saturation adjustment
     989                              qvc(i) = rneb(i,k) * zqs(i)
    859990                          ENDIF
    860991
    861                           rnebss(i,k)=0.0   !--added by OB (needed because of the convergence loop on time)
    862                           fcontrN(i,k)=0.0  !--idem
    863                           fcontrP(i,k)=0.0  !--idem
    864                           qss(i,k)=0.0      !--idem
    865 
    866                        ELSE ! in case of ice supersaturation by Audran
    867 
    868                         !------------------------------------
    869                         ! ICE SUPERSATURATION
    870                         !------------------------------------
    871 
    872                         CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), &
    873                              gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),                 &
    874                              rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                   &
    875                              Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
    876 
    877                         ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))
     992                      ENDIF ! keepgoing
     993                  ENDDO ! loop on klon
     994
     995                  ENDIF ! .NOT. ok_ice_supersat
     996
     997                  DO i=1,klon
     998                      IF (keepgoing(i)) THEN
    878999
    8791000                        ! If vertical heterogeneity, change fraction by volume as well
    … …  
    9581079                    zqn(i) = zq(i)
    9591080                    rneb(i,k) = 1.0
    960                     zcond(i) = MAX(0.0,zqn(i)-zqs(i))
     1081                    IF ( ok_unadjusted_clouds ) THEN
     1082                      !--AB We relax the saturation adjustment assumption
     1083                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
     1084                      zcond(i) = MAX(0., zqn(i) - qvc(i))
     1085                    ELSE
     1086                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
     1087                    ENDIF
    9611088                    rhcl(i,k)=1.0
    9621089                ELSE
    963                     zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
     1090                    IF ( ok_unadjusted_clouds ) THEN
     1091                      !--AB We relax the saturation adjustment assumption
     1092                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
     1093                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
     1094                    ELSE
     1095                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
     1096                    ENDIF
    9641097                    ! following line is very strange and probably wrong:
    9651098                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
    … …  
    9851118          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
    9861119        ENDIF
     1120   
     1121    !--AB Write diagnostics and tracers for ice supersaturation
     1122    IF ( ok_ice_supersat ) THEN
     1123      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
     1124      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
     1125
     1126      DO i = 1, klon
     1127
     1128        cf_seri(i,k) = rneb(i,k)
     1129
     1130        IF ( .NOT. ok_unadjusted_clouds ) THEN
     1131          qvc(i) = zqs(i) * rneb(i,k)
     1132        ENDIF
     1133        IF ( zq(i) .GT. min_qParent ) THEN
     1134          rvc_seri(i,k) = qvc(i) / zq(i)
     1135        ELSE
     1136          rvc_seri(i,k) = min_ratio
     1137        ENDIF
     1138        !--The MIN barrier is NEEDED because of:
     1139        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
     1140        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
     1141        !--The MAX barrier is a safeguard that should not be activated
     1142        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
     1143
     1144        !--Diagnostics
     1145        gamma_cond(i,k) = gammasat(i)
     1146        IF ( issrfra(i,k) .LT. eps ) THEN
     1147          issrfra(i,k) = 0.
     1148          qissr(i,k) = 0.
     1149        ENDIF
     1150        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
     1151        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
     1152        IF ( subfra(i,k) .LT. eps ) THEN
     1153          subfra(i,k) = 0.
     1154          qsub(i,k) = 0.
     1155        ENDIF
     1156        qcld(i,k) = qvc(i) + zcond(i)
     1157        IF ( cf_seri(i,k) .LT. eps ) THEN
     1158          qcld(i,k) = 0.
     1159        ENDIF
     1160      ENDDO
     1161    ENDIF
     1162
    9871163
    9881164    ! ----------------------------------------------------------------
    … …  
    14071583    ENDDO
    14081584
    1409     !--save some variables for ice supersaturation
    1410     !
    1411     DO i = 1, klon
    1412         ! for memory
    1413         rneb_seri(i,k) = rneb(i,k)
    1414 
    1415         ! for diagnostics
    1416         rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
    1417 
    1418         qvc(i,k) = zqs(i) * rneb(i,k)
    1419         qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--added by OB because of pathologic cases with the lognormal
    1420         qcld(i,k) = qvc(i,k) + zcond(i)
    1421    ENDDO
    1422    !q_sat
    1423    CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:))
    1424    CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:))
    1425 
    1426 
    1427 
    14281585    ! Outputs:
    14291586    !-------------------------------

  • LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90

    r4915 r4951  
    66  !--------------------
    77 
    8   REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
    9   !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI)
     8  REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W
     9  !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W)
    1010 
    1111  REAL, SAVE, PROTECTED :: seuil_neb=0.001      ! cloud fraction threshold: a cloud can precipitate when exceeded
    … …  
    136136  !$OMP THREADPRIVATE(dist_liq)
    137137
    138   REAL, SAVE, PROTECTED    :: tresh_cl=0.0                  ! cloud fraction threshold for cloud top search
     138  REAL, SAVE, PROTECTED  :: tresh_cl=0.0                  ! cloud fraction threshold for cloud top search
    139139  !$OMP THREADPRIVATE(tresh_cl)
     140
     141  !--Parameters for condensation and ice supersaturation
     142  LOGICAL, SAVE, PROTECTED :: ok_external_lognormal=.FALSE.  ! if True, the lognormal condensation scheme is calculated in the lmdz_lscp_condensation routine
     143  !$OMP THREADPRIVATE(ok_external_lognormal)
     144
     145  LOGICAL, SAVE, PROTECTED :: ok_ice_supersat=.FALSE.        ! activates the condensation scheme that allows for ice supersaturation
     146  !$OMP THREADPRIVATE(ok_ice_supersat)
     147
     148  LOGICAL, SAVE, PROTECTED :: ok_unadjusted_clouds=.FALSE.   ! if True, relax the saturation adjustment assumption for ice clouds
     149  !$OMP THREADPRIVATE(ok_unadjusted_clouds)
     150
     151  LOGICAL, SAVE, PROTECTED :: ok_weibull_warm_clouds=.FALSE. ! if True, the weibull condensation scheme replaces the lognormal condensation scheme at positive temperatures
     152  !$OMP THREADPRIVATE(ok_weibull_warm_clouds)
     153
     154  INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf=3       ! iflag for the distribution of water inside ice clouds
     155  !$OMP THREADPRIVATE(iflag_cloud_sublim_pdf)
     156
     157  REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp=1./3.            ! [-] shape factor of the gamma distribution of water inside ice clouds
     158  !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
     159 
     160  REAL, SAVE, PROTECTED :: beta_pdf_lscp=2.8e-3              ! []
     161  !$OMP THREADPRIVATE(beta_pdf_lscp)
     162 
     163  REAL, SAVE, PROTECTED :: temp_thresh_pdf_lscp=188.         ! [K] below this temperature, water vapor is homoegeneously distributed in the clear sky region
     164  !$OMP THREADPRIVATE(temp_thresh_pdf_lscp)
     165 
     166  REAL, SAVE, PROTECTED :: rhlmid_pdf_lscp=57.6              ! [%]
     167  !$OMP THREADPRIVATE(rhlmid_pdf_lscp)
     168 
     169  REAL, SAVE, PROTECTED :: k0_pdf_lscp=1.75                  ! [-] minimum value of the shape parameter for the weibull law
     170  !$OMP THREADPRIVATE(k0_pdf_lscp)
     171 
     172  REAL, SAVE, PROTECTED :: kappa_pdf_lscp=0.0147             ! []
     173  !$OMP THREADPRIVATE(kappa_pdf_lscp)
     174 
     175  REAL, SAVE, PROTECTED :: rhl0_pdf_lscp=81.5                ! [%]
     176  !$OMP THREADPRIVATE(rhl0_pdf_lscp)
     177 
     178  REAL, SAVE, PROTECTED :: cond_thresh_pdf_lscp=1.E-7        ! [-]
     179  !$OMP THREADPRIVATE(cond_thresh_pdf_lscp)
     180 
     181  REAL, SAVE, PROTECTED :: a_homofreez=2.349                 ! [-]
     182  !$OMP THREADPRIVATE(a_homofreez)
     183 
     184  REAL, SAVE, PROTECTED :: b_homofreez=259.                  ! [K]
     185  !$OMP THREADPRIVATE(b_homofreez)
     186
     187  REAL, SAVE, PROTECTED :: delta_hetfreez=1.                 ! [-] value between 0 and 1 to simulate for heterogeneous freezing.
     188  !$OMP THREADPRIVATE(delta_hetfreez)
     189 
     190  REAL, SAVE, PROTECTED :: coef_mixing_lscp=1e-7             ! [-]
     191  !$OMP THREADPRIVATE(coef_mixing_lscp)
     192 
     193  REAL, SAVE, PROTECTED :: coef_shear_lscp=0.1               ! [-]
     194  !$OMP THREADPRIVATE(coef_shear_lscp)
     195 
     196  REAL, SAVE, PROTECTED :: chi_mixing_lscp=1.1               ! [-]
     197  !$OMP THREADPRIVATE(chi_mixing_lscp)
     198
     199!  REAL, SAVE, PROTECTED :: contrail_cross_section=200000.
     200!  !$OMP THREADPRIVATE(contrail_cross_section)
     201  !--End of the parameters for condensation and ice supersaturation
    140202
    141203  !--Parameters for poprecip
    … …  
    225287CONTAINS
    226288
    227 SUBROUTINE lscp_ini(dtime,lunout_in,prt_level_in,ok_ice_sursat, iflag_ratqs, fl_cor_ebil_in, RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, &
    228                     RVTMP2_in, RTT_in,RD_in,RG_in,RPI_in)
     289SUBROUTINE lscp_ini(dtime,lunout_in,prt_level_in,ok_ice_supersat_in, iflag_ratqs, fl_cor_ebil_in, &
     290                    RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RVTMP2_in, &
     291                    RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in)
    229292
    230293
    231294   USE ioipsl_getin_p_mod, ONLY : getin_p
    232    USE ice_sursat_mod, ONLY: ice_sursat_init
    233295   USE lmdz_cloudth_ini, ONLY : cloudth_ini
    234296
    235297   REAL, INTENT(IN)      :: dtime
    236298   INTEGER, INTENT(IN)   :: lunout_in,prt_level_in,iflag_ratqs,fl_cor_ebil_in
    237    LOGICAL, INTENT(IN)   :: ok_ice_sursat
     299   LOGICAL, INTENT(IN)   :: ok_ice_supersat_in
    238300
    239301   REAL, INTENT(IN)      :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in
    240    REAL, INTENT(IN)      ::  RVTMP2_in, RTT_in, RD_in, RG_in, RPI_in
     302   REAL, INTENT(IN)      :: RVTMP2_in, RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in
    241303   character (len=20) :: modname='lscp_ini_mod'
    242304   character (len=80) :: abort_message
    … …  
    247309    fl_cor_ebil=fl_cor_ebil_in
    248310
     311    ok_ice_supersat=ok_ice_supersat_in
     312
    249313    RG=RG_in
    250314    RD=RD_in
     315    RV=RV_in
    251316    RCPD=RCPD_in
    252317    RLVTT=RLVTT_in
    … …  
    257322    RVTMP2=RVTMP2_in
    258323    RPI=RPI_in
     324    EPS_W=EPS_W_in
    259325
    260326
    … …  
    297363    CALL getin_p('tresh_cl',tresh_cl)
    298364    CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
     365    ! for poprecip
    299366    CALL getin_p('ok_poprecip',ok_poprecip)
    300367    CALL getin_p('ok_corr_vap_evasub',ok_corr_vap_evasub)
    … …  
    316383    CALL getin_p('snow_fallspeed_clr',snow_fallspeed_clr)
    317384    CALL getin_p('snow_fallspeed_cld',snow_fallspeed_cld)
     385    ! for condensation and ice supersaturation
     386    CALL getin_p('ok_external_lognormal',ok_external_lognormal)
     387    CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds)
     388    CALL getin_p('ok_weibull_warm_clouds',ok_weibull_warm_clouds)
     389    CALL getin_p('iflag_cloud_sublim_pdf',iflag_cloud_sublim_pdf)
     390    CALL getin_p('mu_subl_pdf_lscp',mu_subl_pdf_lscp)
     391    CALL getin_p('beta_pdf_lscp',beta_pdf_lscp)
     392    CALL getin_p('temp_thresh_pdf_lscp',temp_thresh_pdf_lscp)
     393    CALL getin_p('rhlmid_pdf_lscp',rhlmid_pdf_lscp)
     394    CALL getin_p('k0_pdf_lscp',k0_pdf_lscp)
     395    CALL getin_p('kappa_pdf_lscp',kappa_pdf_lscp)
     396    CALL getin_p('rhl0_pdf_lscp',rhl0_pdf_lscp)
     397    CALL getin_p('cond_thresh_pdf_lscp',cond_thresh_pdf_lscp)
     398    CALL getin_p('a_homofreez',a_homofreez)
     399    CALL getin_p('b_homofreez',b_homofreez)
     400    CALL getin_p('delta_hetfreez',delta_hetfreez)
     401    CALL getin_p('coef_mixing_lscp',coef_mixing_lscp)
     402    CALL getin_p('coef_shear_lscp',coef_shear_lscp)
     403    CALL getin_p('chi_mixing_lscp',chi_mixing_lscp)
     404    !CALL getin_p('contrail_cross_section',contrail_cross_section)
    318405
    319406
    … …  
    354441    WRITE(lunout,*) 'lscp_ini, iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp
    355442    WRITE(lunout,*) 'lscp_ini, fl_cor_ebil', fl_cor_ebil
     443    ! for poprecip
    356444    WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
    357445    WRITE(lunout,*) 'lscp_ini, ok_corr_vap_evasub', ok_corr_vap_evasub
    … …  
    367455    WRITE(lunout,*) 'lscp_ini, snow_fallspeed_clr:', snow_fallspeed_clr
    368456    WRITE(lunout,*) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld
     457    ! for condensation and ice supersaturation
     458    WRITE(lunout,*) 'lscp_ini, ok_external_lognormal:', ok_external_lognormal
     459    WRITE(lunout,*) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat
     460    WRITE(lunout,*) 'lscp_ini, ok_unadjusted_clouds:', ok_unadjusted_clouds
     461    WRITE(lunout,*) 'lscp_ini, ok_weibull_warm_clouds:', ok_weibull_warm_clouds
     462    WRITE(lunout,*) 'lscp_ini, iflag_cloud_sublim_pdf:', iflag_cloud_sublim_pdf
     463    WRITE(lunout,*) 'lscp_ini, mu_subl_pdf_lscp:', mu_subl_pdf_lscp
     464    WRITE(lunout,*) 'lscp_ini, beta_pdf_lscp:', beta_pdf_lscp
     465    WRITE(lunout,*) 'lscp_ini, temp_thresh_pdf_lscp:', temp_thresh_pdf_lscp
     466    WRITE(lunout,*) 'lscp_ini, rhlmid_pdf_lscp:', rhlmid_pdf_lscp
     467    WRITE(lunout,*) 'lscp_ini, k0_pdf_lscp:', k0_pdf_lscp
     468    WRITE(lunout,*) 'lscp_ini, kappa_pdf_lscp:', kappa_pdf_lscp
     469    WRITE(lunout,*) 'lscp_ini, rhl0_pdf_lscp:', rhl0_pdf_lscp
     470    WRITE(lunout,*) 'lscp_ini, a_homofreez:', a_homofreez
     471    WRITE(lunout,*) 'lscp_ini, b_homofreez:', b_homofreez
     472    WRITE(lunout,*) 'lscp_ini, delta_hetfreez', delta_hetfreez
     473    WRITE(lunout,*) 'lscp_ini, coef_mixing_lscp:', coef_mixing_lscp
     474    WRITE(lunout,*) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp
     475    WRITE(lunout,*) 'lscp_ini, chi_mixing_lscp:', chi_mixing_lscp
     476!    WRITE(lunout,*) 'lscp_ini, contrail_cross_section:', contrail_cross_section
    369477
    370478
    … …  
    389497
    390498
     499    !--Check flags for condensation and ice supersaturation
     500    IF ( ok_external_lognormal .AND. ok_ice_supersat ) THEN
     501      abort_message = 'in lscp, ok_external_lognormal=y is incompatible with ok_ice_supersat=y'
     502      CALL abort_physic (modname,abort_message,1)
     503    ENDIF
     504
     505    IF ( ok_weibull_warm_clouds .AND. .NOT. ok_ice_supersat ) THEN
     506      abort_message = 'in lscp, ok_weibull_warm_clouds=y needs ok_ice_supersat=y'
     507      CALL abort_physic (modname,abort_message,1)
     508    ENDIF
     509
     510    IF ( ok_unadjusted_clouds .AND. .NOT. ok_ice_supersat ) THEN
     511      abort_message = 'in lscp, ok_unadjusted_clouds=y needs ok_ice_supersat=y'
     512      CALL abort_physic (modname,abort_message,1)
     513    ENDIF
     514
     515
    391516    !AA Temporary initialisation
    392517    a_tr_sca(1) = -0.5
    … …  
    395520    a_tr_sca(4) = -0.5
    396521   
    397     IF (ok_ice_sursat) CALL ice_sursat_init()
    398 
    399522    CALL cloudth_ini(iflag_cloudth_vert,iflag_ratqs)
    400523

  • LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90

    r4818 r4951  
    311311!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    312312
    313     use lmdz_lscp_ini, only: iflag_gammasat, t_glace_min, RTT
     313    use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
     314    use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez
    314315
    315316    IMPLICIT NONE
    … …  
    326327
    327328    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
    328     REAL  fcirrus, fac
    329     REAL, PARAMETER :: acirrus=2.349
    330     REAL, PARAMETER :: bcirrus=259.0
     329    REAL  f_homofreez, fac
    331330
    332331    INTEGER i
    … …  
    342341            dgammasatdt(i)=0.
    343342
    344         ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
     343        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. temp_nowater)) THEN
    345344           
    346345            IF (iflag_gammasat .GE. 2) THEN         
    … …  
    355354
    356355            IF (iflag_gammasat .GE.1) THEN
    357             ! homogeneous freezing of aerosols, according to
    358             ! Koop, 2000 and Karcher 2008, QJRMS
    359             ! 'Cirrus regime'
    360                fcirrus=acirrus-temp(i)/bcirrus
    361                IF (fcirrus .GT. qsl(i)/qsi(i)) THEN
     356               ! homogeneous freezing of aerosols, according to
     357               ! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
     358               ! 'Cirrus regime'
     359               f_homofreez=a_homofreez-temp(i)/b_homofreez
     360               IF (f_homofreez .GT. qsl(i)/qsi(i)) THEN
    362361                  gammasat(i)=qsl(i)/qsi(i)
    363362                  dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
    364363               ELSE
    365                   gammasat(i)=fcirrus
    366                   dgammasatdt(i)=-1.0/bcirrus
     364                  gammasat(i)= (1.-delta_hetfreez) + delta_hetfreez * f_homofreez
     365                  dgammasatdt(i)=-delta_hetfreez/b_homofreez
    367366               ENDIF
    368367           
    … …  
    452451!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    453452
     453!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     454FUNCTION GAMMAINC ( p, x )
     455
     456!*****************************************************************************80
     457!
     458!! GAMMAINC computes the regularized lower incomplete Gamma Integral
     459!
     460!  Modified:
     461!
     462!    20 January 2008
     463!
     464!  Author:
     465!
     466!    Original FORTRAN77 version by B Shea.
     467!    FORTRAN90 version by John Burkardt.
     468!
     469!  Reference:
     470!
     471!    B Shea,
     472!    Algorithm AS 239:
     473!    Chi-squared and Incomplete Gamma Integral,
     474!    Applied Statistics,
     475!    Volume 37, Number 3, 1988, pages 466-473.
     476!
     477!  Parameters:
     478!
     479!    Input, real X, P, the parameters of the incomplete
     480!    gamma ratio.  0 <= X, and 0 < P.
     481!
     482!    Output, real GAMMAINC, the value of the incomplete
     483!    Gamma integral.
     484!
     485  IMPLICIT NONE
     486
     487  REAL A
     488  REAL AN
     489  REAL ARG
     490  REAL B
     491  REAL C
     492  REAL, PARAMETER :: ELIMIT = - 88.0E+00
     493  REAL GAMMAINC
     494  REAL, PARAMETER :: OFLO = 1.0E+37
     495  REAL P
     496  REAL, PARAMETER :: PLIMIT = 1000.0E+00
     497  REAL PN1
     498  REAL PN2
     499  REAL PN3
     500  REAL PN4
     501  REAL PN5
     502  REAL PN6
     503  REAL RN
     504  REAL, PARAMETER :: TOL = 1.0E-14
     505  REAL X
     506  REAL, PARAMETER :: XBIG = 1.0E+08
     507
     508  GAMMAINC = 0.0E+00
     509
     510  IF ( X == 0.0E+00 ) THEN
     511    GAMMAINC = 0.0E+00
     512    RETURN
     513  END IF
     514!
     515!  IF P IS LARGE, USE A NORMAL APPROXIMATION.
     516!
     517  IF ( PLIMIT < P ) THEN
     518
     519    PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
     520    + 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )
     521
     522    GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
     523    RETURN
     524
     525  END IF
     526!
     527!  IF X IS LARGE SET GAMMAD = 1.
     528!
     529  IF ( XBIG < X ) THEN
     530    GAMMAINC = 1.0E+00
     531    RETURN
     532  END IF
     533!
     534!  USE PEARSON'S SERIES EXPANSION.
     535!  (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
     536!
     537  IF ( X <= 1.0E+00 .OR. X < P ) THEN
     538
     539    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
     540    C = 1.0E+00
     541    GAMMAINC = 1.0E+00
     542    A = P
     543
     544    DO
     545
     546      A = A + 1.0E+00
     547      C = C * X / A
     548      GAMMAINC = GAMMAINC + C
     549
     550      IF ( C <= TOL ) THEN
     551        EXIT
     552      END IF
     553
     554    END DO
     555
     556    ARG = ARG + LOG ( GAMMAINC )
     557
     558    IF ( ELIMIT <= ARG ) THEN
     559      GAMMAINC = EXP ( ARG )
     560    ELSE
     561      GAMMAINC = 0.0E+00
     562    END IF
     563!
     564!  USE A CONTINUED FRACTION EXPANSION.
     565!
     566  ELSE
     567
     568    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
     569    A = 1.0E+00 - P
     570    B = A + X + 1.0E+00
     571    C = 0.0E+00
     572    PN1 = 1.0E+00
     573    PN2 = X
     574    PN3 = X + 1.0E+00
     575    PN4 = X * B
     576    GAMMAINC = PN3 / PN4
     577
     578    DO
     579
     580      A = A + 1.0E+00
     581      B = B + 2.0E+00
     582      C = C + 1.0E+00
     583      AN = A * C
     584      PN5 = B * PN3 - AN * PN1
     585      PN6 = B * PN4 - AN * PN2
     586
     587      IF ( PN6 /= 0.0E+00 ) THEN
     588
     589        RN = PN5 / PN6
     590
     591        IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
     592          EXIT
     593        END IF
     594
     595        GAMMAINC = RN
     596
     597      END IF
     598
     599      PN1 = PN3
     600      PN2 = PN4
     601      PN3 = PN5
     602      PN4 = PN6
     603!
     604!  RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
     605!
     606      IF ( OFLO <= ABS ( PN5 ) ) THEN
     607        PN1 = PN1 / OFLO
     608        PN2 = PN2 / OFLO
     609        PN3 = PN3 / OFLO
     610        PN4 = PN4 / OFLO
     611      END IF
     612
     613    END DO
     614
     615    ARG = ARG + LOG ( GAMMAINC )
     616
     617    IF ( ELIMIT <= ARG ) THEN
     618      GAMMAINC = 1.0E+00 - EXP ( ARG )
     619    ELSE
     620      GAMMAINC = 1.0E+00
     621    END IF
     622
     623  END IF
     624
     625  RETURN
     626END FUNCTION GAMMAINC
     627!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     628
    454629END MODULE lmdz_lscp_tools
    455630

  • LMDZ6/branches/cirrus/libf/phylmd/phyetat0_mod.F90

    r4744 r4951  
    2121       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
    2222       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
    23        ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, rneb_ancien, radpas, radsol, rain_fall, ratqs, &
     23       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, &
     24       cf_ancien, rvc_ancien, radpas, radsol, rain_fall, ratqs, &
    2425       rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
    2526       solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
    … …  
    397398  ancien_ok=ancien_ok.AND.phyetat0_get(ql_ancien,"QLANCIEN","QLANCIEN",0.)
    398399  ancien_ok=ancien_ok.AND.phyetat0_get(qs_ancien,"QSANCIEN","QSANCIEN",0.)
    399   ancien_ok=ancien_ok.AND.phyetat0_get(rneb_ancien,"RNEBANCIEN","RNEBANCIEN",0.)
    400400  ancien_ok=ancien_ok.AND.phyetat0_get(u_ancien,"UANCIEN","UANCIEN",0.)
    401401  ancien_ok=ancien_ok.AND.phyetat0_get(v_ancien,"VANCIEN","VANCIEN",0.)
    … …  
    412412     prbsw_ancien(:)=0.
    413413  ENDIF
     414 
     415  ! cas specifique des variables de la sursaturation par rapport a la glace
     416  IF ( ok_ice_supersat ) THEN
     417    ancien_ok=ancien_ok.AND.phyetat0_get(cf_ancien,"CFANCIEN","CFANCIEN",0.)
     418    ancien_ok=ancien_ok.AND.phyetat0_get(rvc_ancien,"RVCANCIEN","RVCANCIEN",0.)
     419  ELSE
     420    cf_ancien(:,:)=0.
     421    rvc_ancien(:,:)=0.
     422  ENDIF
    414423
    415424  ! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
    … …  
    419428       (maxval(ql_ancien).EQ.minval(ql_ancien))     .OR. &
    420429       (maxval(qs_ancien).EQ.minval(qs_ancien))     .OR. &
    421        (maxval(rneb_ancien).EQ.minval(rneb_ancien)) .OR. &
    422430       (maxval(prw_ancien).EQ.minval(prw_ancien))   .OR. &
    423431       (maxval(prlw_ancien).EQ.minval(prlw_ancien)) .OR. &
    … …  
    432440       ancien_ok=.false.
    433441    ENDIF
     442  ENDIF
     443
     444  IF ( ok_ice_supersat ) THEN
     445    IF ( (maxval(cf_ancien).EQ.minval(cf_ancien))     .OR. &
     446         (maxval(rvc_ancien).EQ.minval(rvc_ancien)) ) THEN
     447       ancien_ok=.false.
     448     ENDIF
    434449  ENDIF
    435450

  • LMDZ6/branches/cirrus/libf/phylmd/phyredem.F90

    r4744 r4951  
    1919                                zval, rugoro, t_ancien, q_ancien,            &
    2020                                prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
    21                                 ql_ancien, qs_ancien, qbs_ancien, rneb_ancien, u_ancien, &
    22                                 v_ancien, clwcon, rnebcon, ratqs, pbl_tke,   &
     21                                ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
     22                                rvc_ancien, u_ancien, v_ancien,              &
     23                                clwcon, rnebcon, ratqs, pbl_tke,             &
    2324                                wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
    2425                                wake_deltat, wake_deltaq, wake_s, awake_s,   &
    … …  
    252253    ENDIF
    253254
    254     CALL put_field(pass,"RNEBANCIEN", "RNEBANCIEN", rneb_ancien)
     255    IF ( ok_ice_supersat ) THEN
     256      CALL put_field(pass,"CFANCIEN", "CFANCIEN", cf_ancien)
     257      CALL put_field(pass,"RVCANCIEN", "RVCANCIEN", rvc_ancien)
     258    ENDIF
    255259
    256260    CALL put_field(pass,"PRWANCIEN", "PRWANCIEN", prw_ancien)

  • LMDZ6/branches/cirrus/libf/phylmd/phys_local_var_mod.F90

    r4926 r4951  
    1818      REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
    1919      !$OMP THREADPRIVATE(u_seri, v_seri)
    20       REAL, SAVE, ALLOCATABLE :: rneb_seri(:,:)
    21       !$OMP THREADPRIVATE(rneb_seri)
    22       REAL, SAVE, ALLOCATABLE :: d_rneb_dyn(:,:)
    23       !$OMP THREADPRIVATE(d_rneb_dyn)
     20      REAL, SAVE, ALLOCATABLE :: cf_seri(:,:), rvc_seri(:,:)
     21      !$OMP THREADPRIVATE(cf_seri, rvc_seri)
    2422      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),wprime(:,:,:)
    2523      !$OMP THREADPRIVATE(l_mixmin, l_mix, wprime)
    … …  
    3836      REAL, SAVE, ALLOCATABLE :: d_u_dyn(:,:), d_v_dyn(:,:)
    3937      !$OMP THREADPRIVATE(d_u_dyn, d_v_dyn)
     38      REAL, SAVE, ALLOCATABLE :: d_cf_dyn(:,:), d_rvc_dyn(:,:)
     39      !$OMP THREADPRIVATE(d_cf_dyn, d_rvc_dyn)
    4040      REAL, SAVE, ALLOCATABLE :: d_tr_dyn(:,:,:)
    4141      !$OMP THREADPRIVATE(d_tr_dyn)
    … …  
    494494!$OMP THREADPRIVATE(zn2mout)
    495495
    496       REAL, SAVE, ALLOCATABLE :: qclr(:,:)
    497       !$OMP THREADPRIVATE(qclr)
    498       REAL, SAVE, ALLOCATABLE :: qcld(:,:)
    499       !$OMP THREADPRIVATE(qcld)
    500       REAL, SAVE, ALLOCATABLE :: qss(:,:)
    501       !$OMP THREADPRIVATE(qss)
    502       REAL, SAVE, ALLOCATABLE :: qvc(:,:)
    503       !$OMP THREADPRIVATE(qvc)
    504       REAL, SAVE, ALLOCATABLE :: rnebclr(:,:)
    505       !$OMP THREADPRIVATE(rnebclr)
    506       REAL, SAVE, ALLOCATABLE :: rnebss(:,:)
    507       !$OMP THREADPRIVATE(rnebss)
    508       REAL, SAVE, ALLOCATABLE :: gamma_ss(:,:)
    509       !$OMP THREADPRIVATE(gamma_ss)
    510       REAL, SAVE, ALLOCATABLE :: N1_ss(:,:)
    511       !$OMP THREADPRIVATE(N1_ss)
    512       REAL, SAVE, ALLOCATABLE :: N2_ss(:,:)
    513       !$OMP THREADPRIVATE(N2_ss)
    514       REAL, SAVE, ALLOCATABLE :: drneb_sub(:,:)
    515       !$OMP THREADPRIVATE(drneb_sub)
    516       REAL, SAVE, ALLOCATABLE :: drneb_con(:,:)
    517       !$OMP THREADPRIVATE(drneb_con)
    518       REAL, SAVE, ALLOCATABLE :: drneb_tur(:,:)
    519       !$OMP THREADPRIVATE(drneb_tur)
    520       REAL, SAVE, ALLOCATABLE :: drneb_avi(:,:)
    521       !$OMP THREADPRIVATE(drneb_avi)
    522       REAL, SAVE, ALLOCATABLE :: zqsatl(:,:)
    523       !$OMP THREADPRIVATE(zqsatl)
    524       REAL, SAVE, ALLOCATABLE :: zqsats(:,:)
    525       !$OMP THREADPRIVATE(zqsats)
    526       REAL, SAVE, ALLOCATABLE :: Tcontr(:,:)
    527       !$OMP THREADPRIVATE(Tcontr)
    528       REAL, SAVE, ALLOCATABLE :: qcontr(:,:)
    529       !$OMP THREADPRIVATE(qcontr)
    530       REAL, SAVE, ALLOCATABLE :: qcontr2(:,:)
    531       !$OMP THREADPRIVATE(qcontr2)
    532       REAL, SAVE, ALLOCATABLE :: fcontrN(:,:)
    533       !$OMP THREADPRIVATE(fcontrN)
    534       REAL, SAVE, ALLOCATABLE :: fcontrP(:,:)
    535       !$OMP THREADPRIVATE(fcontrP)
     496!-- LSCP - condensation and ice supersaturation variables
     497      REAL, SAVE, ALLOCATABLE :: qsub(:,:), qissr(:,:), qcld(:,:)
     498      !$OMP THREADPRIVATE(qsub, qissr, qcld)
     499      REAL, SAVE, ALLOCATABLE :: subfra(:,:), issrfra(:,:)
     500      !$OMP THREADPRIVATE(subfra, issrfra)
     501      REAL, SAVE, ALLOCATABLE :: gamma_cond(:,:)
     502      !$OMP THREADPRIVATE(gamma_cond)
     503      REAL, SAVE, ALLOCATABLE :: ratio_qi_qtot(:,:)
     504      !$OMP THREADPRIVATE(ratio_qi_qtot)
     505      REAL, SAVE, ALLOCATABLE :: dcf_sub(:,:), dcf_con(:,:), dcf_mix(:,:)
     506      !$OMP THREADPRIVATE(dcf_sub, dcf_con, dcf_mix)
     507      REAL, SAVE, ALLOCATABLE :: dqi_adj(:,:), dqi_sub(:,:), dqi_con(:,:), dqi_mix(:,:)
     508      !$OMP THREADPRIVATE(dqi_adj, dqi_sub, dqi_con, dqi_mix)
     509      REAL, SAVE, ALLOCATABLE :: dqvc_adj(:,:), dqvc_sub(:,:), dqvc_con(:,:), dqvc_mix(:,:)
     510      !$OMP THREADPRIVATE(dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix)
     511      REAL, SAVE, ALLOCATABLE :: qsatliq(:,:), qsatice(:,:)
     512      !$OMP THREADPRIVATE(qsatliq, qsatice)
     513
     514!-- LSCP - aviation and contrails variables
     515      REAL, SAVE, ALLOCATABLE :: Tcontr(:,:), qcontr(:,:), qcontr2(:,:)
     516      !$OMP THREADPRIVATE(Tcontr, qcontr, qcontr2)
     517      REAL, SAVE, ALLOCATABLE :: fcontrN(:,:), fcontrP(:,:)
     518      !$OMP THREADPRIVATE(fcontrN, fcontrP)
     519      REAL, SAVE, ALLOCATABLE :: dcf_avi(:,:), dqi_avi(:,:), dqvc_avi(:,:)
     520      !$OMP THREADPRIVATE(dcf_avi, dqi_avi, dqvc_avi)
     521      REAL, SAVE, ALLOCATABLE :: flight_dist(:,:), flight_h2o(:,:)
     522      !$OMP THREADPRIVATE(flight_dist, flight_h2o)
     523
     524!-- LSCP - mixed phase clouds variables
    536525      REAL, SAVE, ALLOCATABLE :: distcltop(:,:)
    537526      !$OMP THREADPRIVATE(distcltop)
    … …  
    539528      !$OMP THREADPRIVATE(temp_cltop)
    540529
    541 
    542 !--POPRECIP variables
     530!-- LSCP - POPRECIP variables
    543531      REAL, SAVE, ALLOCATABLE :: qraindiag(:,:)
    544532      !$OMP THREADPRIVATE(qraindiag)
    … …  
    668656      ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev), qbs_seri(klon,klev))
    669657      ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
     658      ALLOCATE(cf_seri(klon,klev),rvc_seri(klon,klev))
    670659      ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
    671660      ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
    … …  
    678667      ALLOCATE(d_q_dyn2d(klon),d_ql_dyn2d(klon),d_qs_dyn2d(klon), d_qbs_dyn2d(klon))
    679668      ALLOCATE(d_u_dyn(klon,klev),d_v_dyn(klon,klev))
     669      ALLOCATE(d_cf_dyn(klon,klev),d_rvc_dyn(klon,klev))
    680670      ALLOCATE(d_tr_dyn(klon,klev,nbtr))                   !RomP
    681671      ALLOCATE(d_t_con(klon,klev),d_q_con(klon,klev),d_q_con_zmasse(klon,klev))
    … …  
    954944      ALLOCATE(zn2mout(klon,6))
    955945
    956 ! Supersaturation
    957       ALLOCATE(rneb_seri(klon,klev))
    958       ALLOCATE(d_rneb_dyn(klon,klev))
    959       ALLOCATE(qclr(klon,klev), qcld(klon,klev), qss(klon,klev), qvc(klon,klev))
    960       ALLOCATE(rnebclr(klon,klev), rnebss(klon,klev), gamma_ss(klon,klev))
    961       ALLOCATE(N1_ss(klon,klev), N2_ss(klon,klev))
    962       ALLOCATE(drneb_sub(klon,klev), drneb_con(klon,klev), drneb_tur(klon,klev), drneb_avi(klon,klev))
    963       ALLOCATE(zqsatl(klon,klev), zqsats(klon,klev))
    964       ALLOCATE(Tcontr(klon,klev), qcontr(klon,klev), qcontr2(klon,klev), fcontrN(klon,klev), fcontrP(klon,klev))
    965 
    966 !--POPRECIP variables
     946!-- LSCP - condensation and ice supersaturation variables
     947      ALLOCATE(qsub(klon,klev), qissr(klon,klev), qcld(klon,klev))
     948      ALLOCATE(subfra(klon,klev), issrfra(klon,klev))
     949      ALLOCATE(gamma_cond(klon,klev), ratio_qi_qtot(klon,klev))
     950      ALLOCATE(dcf_sub(klon,klev), dcf_con(klon,klev), dcf_mix(klon,klev))
     951      ALLOCATE(dqi_adj(klon,klev), dqi_sub(klon,klev), dqi_con(klon,klev), dqi_mix(klon,klev))
     952      ALLOCATE(dqvc_adj(klon,klev), dqvc_sub(klon,klev), dqvc_con(klon,klev), dqvc_mix(klon,klev))
     953      ALLOCATE(qsatliq(klon,klev), qsatice(klon,klev))
     954
     955!-- LSCP - aviation and contrails variables
     956      ALLOCATE(Tcontr(klon,klev), qcontr(klon,klev), qcontr2(klon,klev))
     957      ALLOCATE(fcontrN(klon,klev), fcontrP(klon,klev))
     958      ALLOCATE(dcf_avi(klon,klev), dqi_avi(klon,klev), dqvc_avi(klon,klev))
     959      ALLOCATE(flight_dist(klon,klev), flight_h2o(klon,klev))
     960
     961!-- LSCP - POPRECIP variables
    967962      ALLOCATE(qraindiag(klon,klev), qsnowdiag(klon,klev))
    968963      ALLOCATE(dqreva(klon,klev), dqssub(klon,klev))
    … …  
    10221017      DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri)
    10231018      DEALLOCATE(u_seri,v_seri)
     1019      DEALLOCATE(cf_seri,rvc_seri)
    10241020      DEALLOCATE(l_mixmin,l_mix,wprime)
    10251021      DEALLOCATE(pbl_eps)
    … …  
    10301026      DEALLOCATE(d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, d_qbs_dyn2d)
    10311027      DEALLOCATE(d_u_dyn,d_v_dyn)
     1028      DEALLOCATE(d_cf_dyn,d_rvc_dyn)
    10321029      DEALLOCATE(d_tr_dyn)                      !RomP
    10331030      DEALLOCATE(d_t_con,d_q_con,d_q_con_zmasse)
    … …  
    12701267      DEALLOCATE(zn2mout)
    12711268
    1272 ! Supersaturation
    1273       DEALLOCATE(rneb_seri)
    1274       DEALLOCATE(d_rneb_dyn)
    1275       DEALLOCATE(qclr, qcld, qss, qvc)
    1276       DEALLOCATE(rnebclr, rnebss, gamma_ss)
    1277       DEALLOCATE(N1_ss, N2_ss)
    1278       DEALLOCATE(drneb_sub, drneb_con, drneb_tur, drneb_avi)
    1279       DEALLOCATE(zqsatl, zqsats)
    1280       DEALLOCATE(Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
    1281 
    1282 !--POPRECIP variables
     1269!-- LSCP - condensation and ice supersaturation variables
     1270      DEALLOCATE(qsub, qissr, qcld)
     1271      DEALLOCATE(subfra, issrfra)
     1272      DEALLOCATE(gamma_cond, ratio_qi_qtot)
     1273      DEALLOCATE(dcf_sub, dcf_con, dcf_mix)
     1274      DEALLOCATE(dqi_adj, dqi_sub, dqi_con, dqi_mix)
     1275      DEALLOCATE(dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix)
     1276      DEALLOCATE(qsatliq, qsatice)
     1277
     1278!-- LSCP - aviation and contrails variables
     1279      DEALLOCATE(Tcontr, qcontr, qcontr2)
     1280      DEALLOCATE(fcontrN, fcontrP)
     1281      DEALLOCATE(dcf_avi, dqi_avi, dqvc_avi)
     1282      DEALLOCATE(flight_dist, flight_h2o)
     1283
     1284!-- LSCP - POPRECIP variables
    12831285      DEALLOCATE(qraindiag, qsnowdiag)
    12841286      DEALLOCATE(dqreva, dqssub)

  • LMDZ6/branches/cirrus/libf/phylmd/phys_output_ctrlout_mod.F90

    r4887 r4951  
    450450  TYPE(ctrl_out), SAVE :: o_SWdn200clr = ctrl_out((/ 10, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
    451451    'SWdn200clr', 'SWdn clear sky at 200mb', 'W/m2', (/ ('', i=1, 10) /))
    452 
    453   ! arajouter
    454   !  type(ctrl_out),save :: o_LWupTOA     = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'LWupTOA', &
    455   !    (/ ('', i=1, 10) /))
    456   !  type(ctrl_out),save :: o_LWupTOAclr  = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'LWupTOAclr', &
    457   !    (/ ('', i=1, 10) /))
     452  TYPE(ctrl_out), SAVE :: o_LWupTOA = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/), &
     453    'LWupTOA', 'LWup at TOA', 'W/m2', (/ ('', i=1, 10) /))
     454  TYPE(ctrl_out), SAVE :: o_LWupTOAclr = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/), &
     455    'LWupTOAclr', 'LWup clear sky at TOA', 'W/m2', (/ ('', i=1, 10) /))
    458456  !  type(ctrl_out),save :: o_LWdnTOA     = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'LWdnTOA', &
    459457  !    (/ ('', i=1, 10) /))
    … …  
    20692067    'albslw3', 'Surface albedo LW3', '-', (/ ('', i=1, 10) /))
    20702068
    2071 !--aviation & supersaturation
    2072   TYPE(ctrl_out), SAVE :: o_oclr = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2073     'oclr', 'Clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
    2074   TYPE(ctrl_out), SAVE :: o_ocld = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2075     'ocld', 'Cloudy sky total water', 'kg/kg', (/ ('', i=1, 10) /))
    2076   TYPE(ctrl_out), SAVE :: o_oss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2077     'oss', 'ISSR total water', 'kg/kg', (/ ('', i=1, 10) /))
    2078   TYPE(ctrl_out), SAVE :: o_ovc = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2079     'ovc', 'In-cloup vapor', 'kg/kg', (/ ('', i=1, 10) /))
    2080   TYPE(ctrl_out), SAVE :: o_rnebclr = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2081     'rnebclr', 'Clear sky fraction', '-', (/ ('', i=1, 10) /))
    2082   TYPE(ctrl_out), SAVE :: o_rnebss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2083     'rnebss', 'ISSR fraction', '-', (/ ('', i=1, 10) /))
    2084   TYPE(ctrl_out), SAVE :: o_rnebseri = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2085     'rnebseri', 'Cloud fraction', '-', (/ ('', i=1, 10) /))
    2086   TYPE(ctrl_out), SAVE :: o_gammass = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2087     'gammass', 'Gamma supersaturation', '', (/ ('', i=1, 10) /))
    2088   TYPE(ctrl_out), SAVE :: o_N1_ss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2089     'N1ss', 'N1', '', (/ ('', i=1, 10) /))
    2090   TYPE(ctrl_out), SAVE :: o_N2_ss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2091     'N2ss', 'N2', '', (/ ('', i=1, 10) /))
    2092   TYPE(ctrl_out), SAVE :: o_drnebsub = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2093     'drnebsub', 'Cloud fraction change because of sublimation', 's-1', (/ ('', i=1, 10) /))
    2094   TYPE(ctrl_out), SAVE :: o_drnebcon = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2095     'drnebcon', 'Cloud fraction change because of condensation', 's-1', (/ ('', i=1, 10) /))
    2096   TYPE(ctrl_out), SAVE :: o_drnebtur = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2097     'drnebtur', 'Cloud fraction change because of turbulence', 's-1', (/ ('', i=1, 10) /))
    2098   TYPE(ctrl_out), SAVE :: o_drnebavi = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2099     'drnebavi', 'Cloud fraction change because of aviation', 's-1', (/ ('', i=1, 10) /))
    2100   TYPE(ctrl_out), SAVE :: o_qsatl = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2101     'qsatl', 'Saturation with respect to liquid water', '', (/ ('', i=1, 10) /))
    2102   TYPE(ctrl_out), SAVE :: o_qsats = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2103     'qsats', 'Saturation with respect to solid water', '', (/ ('', i=1, 10) /))
    2104   TYPE(ctrl_out), SAVE :: o_flight_m = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2105     'flightm', 'Flown meters', 'm/s/mesh', (/ ('', i=1, 10) /))
    2106   TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
    2107     'flighth2o', 'H2O flight emission', 'kg H2O/s/mesh', (/ ('', i=1, 10) /))
    2108   TYPE(ctrl_out), SAVE :: o_Tcontr = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
     2069!-- LSCP - condensation and ice supersaturation variables
     2070  TYPE(ctrl_out), SAVE :: o_cfseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2071    'cfseri', 'Cloud fraction', '-', (/ ('', i=1, 10) /))
     2072  TYPE(ctrl_out), SAVE :: o_dcfdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2073    'dcfdyn', 'Dynamics cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2074  TYPE(ctrl_out), SAVE :: o_rvcseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2075    'rvcseri', 'Cloudy water vapor to total water vapor ratio', '-', (/ ('', i=1, 10) /))
     2076  TYPE(ctrl_out), SAVE :: o_drvcdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2077    'drvcdyn', 'Dynamics cloudy water vapor to total water vapor ratio tendency', 's-1', (/ ('', i=1, 10) /))
     2078  TYPE(ctrl_out), SAVE :: o_qsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2079    'qsub', 'Subsaturated clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
     2080  TYPE(ctrl_out), SAVE :: o_qissr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2081    'qissr', 'Supersaturated clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
     2082  TYPE(ctrl_out), SAVE :: o_qcld = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2083    'qcld', 'Cloudy sky total water', 'kg/kg', (/ ('', i=1, 10) /))
     2084  TYPE(ctrl_out), SAVE :: o_subfra = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2085    'subfra', 'Subsaturated clear sky fraction', '-', (/ ('', i=1, 10) /))
     2086  TYPE(ctrl_out), SAVE :: o_issrfra = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2087    'issrfra', 'Supersaturated clear sky fraction', '-', (/ ('', i=1, 10) /))
     2088  TYPE(ctrl_out), SAVE :: o_gammacond = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2089    'gammacond', 'Condensation threshold w.r.t. saturation', '-', (/ ('', i=1, 10) /))
     2090  TYPE(ctrl_out), SAVE :: o_dcfsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2091    'dcfsub', 'Sublimation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2092  TYPE(ctrl_out), SAVE :: o_dcfcon = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2093    'dcfcon', 'Condensation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2094  TYPE(ctrl_out), SAVE :: o_dcfmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2095    'dcfmix', 'Cloud mixing cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2096  TYPE(ctrl_out), SAVE :: o_dqiadj = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2097    'dqiadj', 'Temperature adjustment ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2098  TYPE(ctrl_out), SAVE :: o_dqisub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2099    'dqisub', 'Sublimation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2100  TYPE(ctrl_out), SAVE :: o_dqicon = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2101    'dqicon', 'Condensation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2102  TYPE(ctrl_out), SAVE :: o_dqimix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2103    'dqimix', 'Cloud mixing ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2104  TYPE(ctrl_out), SAVE :: o_dqvcadj = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2105    'dqvcadj', 'Temperature adjustment cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2106  TYPE(ctrl_out), SAVE :: o_dqvcsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2107    'dqvcsub', 'Sublimation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2108  TYPE(ctrl_out), SAVE :: o_dqvccon = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2109    'dqvccon', 'Condensation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2110  TYPE(ctrl_out), SAVE :: o_dqvcmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2111    'dqvcmix', 'Cloud mixing cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2112  TYPE(ctrl_out), SAVE :: o_qsatl = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2113    'qsatl', 'Saturation with respect to liquid', 'kg/kg', (/ ('', i=1, 10) /))
     2114  TYPE(ctrl_out), SAVE :: o_qsati = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2115    'qsati', 'Saturation with respect to ice', 'kg/kg', (/ ('', i=1, 10) /))
     2116
     2117!-- LSCP - aviation variables
     2118  TYPE(ctrl_out), SAVE :: o_Tcontr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    21092119    'Tcontr', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))
    2110   TYPE(ctrl_out), SAVE :: o_qcontr = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
    2111     'qcontr', 'Specific humidity threshold for contrail formation','Pa', (/ ('', i=1, 10) /))
    2112   TYPE(ctrl_out), SAVE :: o_qcontr2 = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
    2113     'qcontr2', 'Specific humidity threshold for contrail formation','kg/kg', (/ ('', i=1, 10) /))
    2114   TYPE(ctrl_out), SAVE :: o_fcontrN = ctrl_out((/ 2, 2, 2, 2, 2, 2, 11, 11, 11, 11/),&
     2120  TYPE(ctrl_out), SAVE :: o_qcontr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2121    'qcontr', 'Specific humidity threshold for contrail formation', 'Pa', (/ ('', i=1, 10) /))
     2122  TYPE(ctrl_out), SAVE :: o_qcontr2 = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2123    'qcontr2', 'Specific humidity threshold for contrail formation', 'kg/kg', (/ ('', i=1, 10) /))
     2124  TYPE(ctrl_out), SAVE :: o_fcontrN = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    21152125    'fcontrN', 'Fraction with non-persistent contrail in clear-sky', '-', (/ ('', i=1,10)/))
    2116   TYPE(ctrl_out), SAVE :: o_fcontrP = ctrl_out((/ 2, 2, 2, 2, 2, 2, 11, 11, 11, 11/),&
     2126  TYPE(ctrl_out), SAVE :: o_fcontrP = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    21172127    'fcontrP', 'Fraction with persistent contrail in ISSR', '-', (/ ('', i=1,10)/))
     2128  TYPE(ctrl_out), SAVE :: o_dcfavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2129    'dcfavi', 'Aviation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2130  TYPE(ctrl_out), SAVE :: o_dqiavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2131    'dqiavi', 'Aviation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2132  TYPE(ctrl_out), SAVE :: o_dqvcavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2133    'dqvcavi', 'Aviation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2134  TYPE(ctrl_out), SAVE :: o_flight_dist = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2135    'flightdist', 'Aviation flown distance', 'm/s/mesh', (/ ('', i=1, 10) /))
     2136  TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2137    'flighth2o', 'Aviation H2O flight emission', 'kg H2O/s/mesh', (/ ('', i=1, 10) /))
    21182138
    21192139!!!!!!!!!!!!! Sorties niveaux standards de pression NMC

  • LMDZ6/branches/cirrus/libf/phylmd/phys_output_write_mod.F90

    r4887 r4951  
    5454         o_SWdnTOAclr, o_nettop, o_SWup200, &
    5555         o_SWup200clr, o_SWdn200, o_SWdn200clr, &
     56         o_LWupTOA, o_LWupTOAclr, &
    5657         o_LWup200, o_LWup200clr, o_LWdn200, &
    5758         o_LWdn200clr, o_sols, o_sols0, &
    … …  
    216217         o_p_tropopause, o_z_tropopause, o_t_tropopause,  &
    217218         o_col_O3_strato, o_col_O3_tropo,                 &
    218 !--aviation & supersaturation
    219          o_oclr, o_ocld, o_oss, o_ovc, o_rnebss, o_rnebclr, o_rnebseri, o_gammass, &
    220          o_N1_ss, o_N2_ss, o_qsatl, o_qsats, &
    221          o_drnebsub, o_drnebcon, o_drnebtur, o_drnebavi, o_flight_m, o_flight_h2o, &
     219!-- LSCP - condensation and ice supersaturation variables
     220         o_cfseri, o_dcfdyn, o_rvcseri, o_drvcdyn, &
     221         o_qsub, o_qissr, o_qcld, o_subfra, o_issrfra, o_gammacond, &
     222         o_dcfsub, o_dcfcon, o_dcfmix, o_dqiadj, o_dqisub, o_dqicon, o_dqimix, &
     223         o_dqvcadj, o_dqvcsub, o_dqvccon, o_dqvcmix, o_qsatl, o_qsati, &
     224!-- LSCP - aviation variables
    222225         o_Tcontr, o_qcontr, o_qcontr2, o_fcontrN, o_fcontrP, &
     226         o_dcfavi, o_dqiavi, o_dqvcavi, o_flight_dist, o_flight_h2o, &
    223227!--interactive CO2
    224228         o_flx_co2_ocean, o_flx_co2_ocean_cor, &
    … …  
    258262#endif
    259263
    260     USE ice_sursat_mod, ONLY: flight_m, flight_h2o
    261264    USE lmdz_lscp_ini, ONLY: ok_poprecip
    262265
    … …  
    323326         cldq, flwp, fiwp, ue, ve, uq, vq, &
    324327         uwat, vwat, &
    325          rneb_seri, d_rneb_dyn, &
    326328         plcl, plfc, wbeff, convoccur, upwd, dnwd, dnwd0, prw, prlw, prsw, prbsw, water_budget, &
    327329         s_pblh, s_pblt, s_lcl, s_therm, uwriteSTD, &
    … …  
    337339         wdtrainA, wdtrainS, wdtrainM, n2, s2, strig, zcong, zlcl_th, proba_notrig, &
    338340         random_notrig, &
    339          qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
    340          N1_ss, N2_ss, zqsatl, zqsats, &
     341         cf_seri, d_cf_dyn, rvc_seri, d_rvc_dyn, &
     342         qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
     343         dcf_sub, dcf_con, dcf_mix, &
     344         dqi_adj, dqi_sub, dqi_con, dqi_mix, &
     345         dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
     346         qsatliq, qsatice, &
    341347         Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
    342          drneb_sub, drneb_con, drneb_tur, drneb_avi, &
     348         dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
    343349         alp_bl_det, alp_bl_fluct_m, alp_bl_conv, &
    344350         alp_bl_stat, alp_bl_fluct_tke, slab_wfbils, &
    … …  
    10661072       CALL histwrite_phy(o_LWupSFCclr, zx_tmp_fi2d)
    10671073       CALL histwrite_phy(o_LWdnSFCclr, sollwdownclr)
     1074       
     1075       IF (vars_defined) THEN
     1076          zx_tmp_fi2d(:) = lwup(:,klevp1)
     1077       ENDIF
     1078       CALL histwrite_phy(o_LWupTOA, zx_tmp_fi2d)
     1079       
     1080       IF (vars_defined) THEN
     1081          zx_tmp_fi2d(:) = lwup0(:,klevp1)
     1082       ENDIF
     1083       CALL histwrite_phy(o_LWupTOAclr, zx_tmp_fi2d)
    10681084
    10691085       IF (type_trac/='inca' .OR. config_inca=='aeNP') THEN
    … …  
    20062022       ENDIF
    20072023
    2008 !--aviation & supersaturation
    2009        IF (ok_ice_sursat) THEN
    2010          CALL histwrite_phy(o_oclr, qclr)
    2011          CALL histwrite_phy(o_ocld, qcld)
    2012          CALL histwrite_phy(o_oss, qss)
    2013          CALL histwrite_phy(o_ovc, qvc)
    2014          CALL histwrite_phy(o_rnebclr, rnebclr)
    2015          CALL histwrite_phy(o_rnebss, rnebss)
    2016          CALL histwrite_phy(o_rnebseri, rneb_seri)
    2017          CALL histwrite_phy(o_gammass, gamma_ss)
    2018          CALL histwrite_phy(o_N1_ss, N1_ss)
    2019          CALL histwrite_phy(o_N2_ss, N2_ss)
    2020          CALL histwrite_phy(o_drnebsub, drneb_sub)
    2021          CALL histwrite_phy(o_drnebcon, drneb_con)
    2022          CALL histwrite_phy(o_drnebtur, drneb_tur)
    2023          CALL histwrite_phy(o_drnebavi, drneb_avi)
    2024          CALL histwrite_phy(o_qsatl, zqsatl)
    2025          CALL histwrite_phy(o_qsats, zqsats)
     2024!-- LSCP - condensation and supersaturation variables
     2025       IF (ok_ice_supersat) THEN
     2026         CALL histwrite_phy(o_cfseri, cf_seri)
     2027         CALL histwrite_phy(o_dcfdyn, d_cf_dyn)
     2028         CALL histwrite_phy(o_rvcseri, rvc_seri)
     2029         CALL histwrite_phy(o_drvcdyn, d_rvc_dyn)
     2030         CALL histwrite_phy(o_qsub, qsub)
     2031         CALL histwrite_phy(o_qissr, qissr)
     2032         CALL histwrite_phy(o_qcld, qcld)
     2033         CALL histwrite_phy(o_subfra, subfra)
     2034         CALL histwrite_phy(o_issrfra, issrfra)
     2035         CALL histwrite_phy(o_gammacond, gamma_cond)
     2036         CALL histwrite_phy(o_dcfsub, dcf_sub)
     2037         CALL histwrite_phy(o_dcfcon, dcf_con)
     2038         CALL histwrite_phy(o_dcfmix, dcf_mix)
     2039         CALL histwrite_phy(o_dqiadj, dqi_adj)
     2040         CALL histwrite_phy(o_dqisub, dqi_sub)
     2041         CALL histwrite_phy(o_dqicon, dqi_con)
     2042         CALL histwrite_phy(o_dqimix, dqi_mix)
     2043         CALL histwrite_phy(o_dqvcadj, dqvc_adj)
     2044         CALL histwrite_phy(o_dqvcsub, dqvc_sub)
     2045         CALL histwrite_phy(o_dqvccon, dqvc_con)
     2046         CALL histwrite_phy(o_dqvcmix, dqvc_mix)
     2047         CALL histwrite_phy(o_qsatl, qsatliq)
     2048         CALL histwrite_phy(o_qsati, qsatice)
     2049       ENDIF
     2050!-- LSCP - aviation variables
     2051       IF (ok_plane_contrail) THEN
    20262052         CALL histwrite_phy(o_Tcontr, Tcontr)
    20272053         CALL histwrite_phy(o_qcontr, qcontr)
    … …  
    20292055         CALL histwrite_phy(o_fcontrN, fcontrN)
    20302056         CALL histwrite_phy(o_fcontrP, fcontrP)
    2031        ENDIF
    2032        IF (ok_plane_contrail) THEN
    2033          CALL histwrite_phy(o_flight_m, flight_m)
     2057         CALL histwrite_phy(o_dcfavi, dcf_avi)
     2058         CALL histwrite_phy(o_dqiavi, dqi_avi)
     2059         CALL histwrite_phy(o_dqvcavi, dqvc_avi)
     2060         CALL histwrite_phy(o_flight_dist, flight_dist)
    20342061       ENDIF
    20352062       IF (ok_plane_h2o) THEN

  • LMDZ6/branches/cirrus/libf/phylmd/phys_state_var_mod.F90

    r4933 r4951  
    9292      REAL, ALLOCATABLE, SAVE :: u_ancien(:,:), v_ancien(:,:)
    9393!$OMP THREADPRIVATE(u_ancien, v_ancien)
     94      REAL, ALLOCATABLE, SAVE :: cf_ancien(:,:), rvc_ancien(:,:)
     95!$OMP THREADPRIVATE(cf_ancien, rvc_ancien)
    9496!!! RomP >>>
    9597      REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
    … …  
    100102      REAL, ALLOCATABLE, SAVE :: clwcon(:,:),rnebcon(:,:)
    101103!$OMP THREADPRIVATE(clwcon,rnebcon)
    102       REAL, ALLOCATABLE, SAVE :: rneb_ancien(:,:)
    103 !$OMP THREADPRIVATE(rneb_ancien)
    104104      REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:),detrain_cv(:,:),fm_cv(:,:)
    105105!$OMP THREADPRIVATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
    … …  
    586586      ALLOCATE(prw_ancien(klon), prlw_ancien(klon), prsw_ancien(klon), prbsw_ancien(klon))
    587587      ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
     588      ALLOCATE(cf_ancien(klon,klev), rvc_ancien(klon,klev))
    588589!!! Rom P >>>
    589590      ALLOCATE(tr_ancien(klon,klev,nbtr))
    590591!!! Rom P <<<
    591592      ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
    592       ALLOCATE(rneb_ancien(klon,klev))
    593593      ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev),detrain_cv(klon,klev),fm_cv(klon,klev))
    594594      ALLOCATE(ratqs(klon,klev))
    … …  
    809809      DEALLOCATE(zthe, zpic, zval)
    810810      DEALLOCATE(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
    811       DEALLOCATE(qs_ancien, ql_ancien, qbs_ancien, rneb_ancien)
     811      DEALLOCATE(qs_ancien, ql_ancien, qbs_ancien)
    812812      DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien)
    813813      DEALLOCATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
    814814      DEALLOCATE(u_ancien, v_ancien)
     815      DEALLOCATE(cf_ancien, rvc_ancien)
    815816      DEALLOCATE(tr_ancien)                           !RomP
    816817      DEALLOCATE(ratqs, pbl_tke,coefh,coefm)

  • LMDZ6/branches/cirrus/libf/phylmd/physiq_mod.F90

    r4887 r4951  
    7373    USE tracinca_mod, ONLY: config_inca
    7474    USE tropopause_m,     ONLY: dyn_tropopause
    75     USE ice_sursat_mod,  ONLY: flight_init, airplane
    7675    USE vampir
    7776    USE write_field_phy
    … …  
    157156       ! [Variables internes non sauvegardees de la physique]
    158157       ! Variables locales pour effectuer les appels en serie
    159        t_seri,q_seri,ql_seri,qs_seri,qbs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
     158       t_seri,q_seri,ql_seri,qs_seri,qbs_seri, &
     159       u_seri,v_seri,cf_seri,rvc_seri,tr_seri, &
    160160       rhcl, &       
    161161       ! Dynamic tendencies (diagnostics)
    162        d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
     162       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn, &
     163       d_u_dyn,d_v_dyn,d_cf_dyn,d_rvc_dyn,d_tr_dyn, &
    163164       d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d,d_qbs_dyn2d, &
    164165       ! Physic tendencies
    … …  
    334335       pfraclr,pfracld, &
    335336       distcltop,temp_cltop, &
    336        zqsatl, zqsats, &
    337        qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
     337       !-- LSCP - condensation and ice supersaturation variables
     338       qsub, qissr, qcld, subfra, issrfra, gamma_cond, ratio_qi_qtot, &
     339       dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
     340       dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
     341       !-- LSCP - aviation and contrails variables
    338342       Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
     343       dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
     344       !
    339345       cldemi,  &
    340346       cldfra, cldtau, fiwc,  &
    … …  
    511517    !
    512518    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
    513     INTEGER,SAVE :: ivap, iliq, isol, irneb, ibs
    514 !$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
     519    INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, irvc
     520!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, irvc)
    515521    !
    516522    !
    … …  
    13621368       iliq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
    13631369       isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
    1364        irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
    13651370       ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
     1371       icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
     1372       irvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
    13661373!       CALL init_etat0_limit_unstruct
    13671374!       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
    … …  
    14161423       ENDIF
    14171424
    1418        IF (ok_ice_sursat.AND.(iflag_ice_thermo.EQ.0)) THEN
    1419           WRITE (lunout, *) ' ok_ice_sursat=y requires iflag_ice_thermo=1 as well'
     1425       IF (ok_ice_supersat.AND.(iflag_ice_thermo.EQ.0)) THEN
     1426          WRITE (lunout, *) ' ok_ice_supersat=y requires iflag_ice_thermo=1 as well'
    14201427          abort_message='see above'
    14211428          CALL abort_physic(modname,abort_message,1)
    14221429       ENDIF
    14231430
    1424        IF (ok_ice_sursat.AND.(nqo.LT.4)) THEN
    1425           WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
    1426                '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
     1431       IF (ok_ice_supersat.AND.(nqo.LT.5)) THEN
     1432          WRITE (lunout, *) ' ok_ice_supersat=y requires 5 H2O tracers ', &
     1433               '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c) but nqo=', nqo, '. Might as well stop here.'
    14271434          abort_message='see above'
    14281435          CALL abort_physic(modname,abort_message,1)
    14291436       ENDIF
    14301437
    1431        IF (ok_plane_h2o.AND..NOT.ok_ice_sursat) THEN
    1432           WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_sursat=y '
     1438       IF (ok_plane_h2o.AND..NOT.ok_ice_supersat) THEN
     1439          WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_supersat=y '
    14331440          abort_message='see above'
    14341441          CALL abort_physic(modname,abort_message,1)
    14351442       ENDIF
    14361443
    1437        IF (ok_plane_contrail.AND..NOT.ok_ice_sursat) THEN
    1438           WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_sursat=y '
     1444       IF (ok_plane_contrail.AND..NOT.ok_ice_supersat) THEN
     1445          WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_supersat=y '
    14391446          abort_message='see above'
    14401447          CALL abort_physic(modname,abort_message,1)
    … …  
    14421449
    14431450        IF (ok_bs) THEN
    1444          IF ((ok_ice_sursat.AND.nqo .LT.5).OR.(.NOT.ok_ice_sursat.AND.nqo.LT.4)) THEN
     1451         IF ((ok_ice_supersat.AND.nqo .LT.6).OR.(.NOT.ok_ice_supersat.AND.nqo.LT.4)) THEN
    14451452             WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
    14461453                               'but nqo=', nqo
    … …  
    18701877   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
    18711878       CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT)
    1872        CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI)
     1879       CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_supersat,iflag_ratqs,fl_cor_ebil, &
     1880                     RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RV,RG,RPI,EPS_W)
    18731881       CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
    18741882                             RVTMP2, RTT,RD,RG, RV, RPI)
    … …  
    23422350                    sollwdown(:))
    23432351
     2352      !--Init for LSCP - condensation
     2353      ratio_qi_qtot(:,:) = 0.
    23442354
    23452355
    … …  
    24392449          q_seri(i,k)  = qx(i,k,ivap)
    24402450          ql_seri(i,k) = qx(i,k,iliq)
    2441           qbs_seri(i,k) = 0.
     2451          qbs_seri(i,k)= 0.
     2452          cf_seri(i,k) = 0.
     2453          rvc_seri(i,k)= 0.
    24422454          !CR: ATTENTION, on rajoute la variable glace
    24432455          IF (nqo.EQ.2) THEN             !--vapour and liquid only
    24442456             qs_seri(i,k) = 0.
    2445              rneb_seri(i,k) = 0.
    24462457          ELSE IF (nqo.EQ.3) THEN        !--vapour, liquid and ice
    24472458             qs_seri(i,k) = qx(i,k,isol)
    2448              rneb_seri(i,k) = 0.
    2449           ELSE IF (nqo.GE.4) THEN        !--vapour, liquid, ice and rneb and blowing snow
     2459          ELSE IF (nqo.GE.4) THEN        !--vapour, liquid, ice, blowing snow, cloud fraction and cloudy water vapor to total water vapor ratio
    24502460             qs_seri(i,k) = qx(i,k,isol)
    2451              IF (ok_ice_sursat) THEN
    2452                rneb_seri(i,k) = qx(i,k,irneb)
     2461             IF (ok_ice_supersat) THEN
     2462               cf_seri(i,k) = qx(i,k,icf)
     2463               rvc_seri(i,k) = qx(i,k,irvc)
    24532464             ENDIF
    24542465             IF (ok_bs) THEN
    … …  
    25262537       d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
    25272538       d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
    2528        d_qbs_dyn(:,:) = (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
     2539       d_qbs_dyn(:,:)= (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
     2540       d_cf_dyn(:,:) = (cf_seri(:,:)-cf_ancien(:,:))/phys_tstep
     2541       d_rvc_dyn(:,:)= (rvc_seri(:,:)-rvc_ancien(:,:))/phys_tstep
    25292542       CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
    25302543       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
    … …  
    25382551       IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
    25392552       ! !! RomP <<<
    2540        !!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
    2541        d_rneb_dyn(:,:)=0.0
    25422553    ELSE
    25432554       d_u_dyn(:,:)  = 0.0
    … …  
    25472558       d_ql_dyn(:,:) = 0.0
    25482559       d_qs_dyn(:,:) = 0.0
     2560       d_qbs_dyn(:,:)= 0.0
     2561       d_cf_dyn(:,:) = 0.0
     2562       d_rvc_dyn(:,:)= 0.0
    25492563       d_q_dyn2d(:)  = 0.0
    25502564       d_ql_dyn2d(:) = 0.0
    … …  
    25542568       IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
    25552569       ! !! RomP <<<
    2556        d_rneb_dyn(:,:)=0.0
    2557        d_qbs_dyn(:,:)=0.0
    25582570       ancien_ok = .TRUE.
    25592571    ENDIF
    … …  
    26642676       ENDIF
    26652677    ENDIF
     2678
     2679    !-- Needed for LSCP - condensation and ice supersaturation
     2680    IF (ok_ice_supersat) THEN
     2681      DO k = 1, klev
     2682        DO i = 1, klon
     2683          IF ( ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) ) .GT. 0. ) THEN
     2684            ratio_qi_qtot(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2685            rvc_seri(i,k) = rvc_seri(i,k) * q_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2686          ELSE
     2687            ratio_qi_qtot(i,k) = 0.
     2688            rvc_seri(i,k) = 0.
     2689          ENDIF
     2690        ENDDO
     2691      ENDDO
     2692    ENDIF
     2693
    26662694    !
    26672695    ! Re-evaporer l'eau liquide nuageuse
    … …  
    38883916
    38893917    !--mise à jour de flight_m et flight_h2o dans leur module
    3890     IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
    3891       CALL airplane(debut,pphis,pplay,paprs,t_seri)
    3892     ENDIF
     3918    !IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
     3919    !  CALL airplane(debut,pphis,pplay,paprs,t_seri)
     3920    !ENDIF
    38933921
    38943922    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
    38953923         t_seri, q_seri,ptconv,ratqs, &
    3896          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
     3924         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
    38973925         pfraclr,pfracld, &
    38983926         radocond, picefra, rain_lsc, snow_lsc, &
    … …  
    39003928         prfl, psfl, rhcl,  &
    39013929         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
    3902          iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, temp_cltop, &
    3903          qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
    3904          Tcontr, qcontr, qcontr2, fcontrN, fcontrP , &
     3930         iflag_ice_thermo, distcltop, temp_cltop, cell_area, &
     3931         cf_seri, rvc_seri, u_seri, v_seri, pbl_eps(:,:,is_ave), &
     3932         qsub, qissr, qcld, subfra, issrfra, gamma_cond, ratio_qi_qtot, &
     3933         dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
     3934         dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
     3935         Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
     3936         dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
    39053937         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
    39063938         qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
    … …  
    56235655             d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
    56245656          ENDIF
    5625           !--ice_sursat: nqo=4, on ajoute rneb
    5626           IF (nqo.ge.4 .and. ok_ice_sursat) THEN
    5627              d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
     5657          !--ice_supersat: nqo=5, we add cloud fraction and cloudy water vapor to total water vapor ratio
     5658          IF (nqo.ge.5 .and. ok_ice_supersat) THEN
     5659             d_qx(i,k,icf) = ( cf_seri(i,k) - qx(i,k,icf) ) / phys_tstep
     5660             d_qx(i,k,irvc) = ( rvc_seri(i,k) - qx(i,k,irvc) ) / phys_tstep
    56285661          ENDIF
    56295662
    … …  
    56595692    ql_ancien(:,:) = ql_seri(:,:)
    56605693    qs_ancien(:,:) = qs_seri(:,:)
    5661     qbs_ancien(:,:) = qbs_seri(:,:)
    5662     rneb_ancien(:,:) = rneb_seri(:,:)
     5694    qbs_ancien(:,:)= qbs_seri(:,:)
     5695    cf_ancien(:,:) = cf_seri(:,:)
     5696    rvc_ancien(:,:)= rvc_seri(:,:)
    56635697    CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
    56645698    CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
Note: See TracChangeset for help on using the changeset viewer.