Index: LMDZ6/branches/cirrus/DefLists/field_def_lmdz.xml
===================================================================
--- LMDZ6/branches/cirrus/DefLists/field_def_lmdz.xml (revision 4950)
+++ LMDZ6/branches/cirrus/DefLists/field_def_lmdz.xml (revision 4951)
@@ -156,4 +156,6 @@
+
+
@@ -798,27 +800,37 @@
pres*ovap*461.5 / (287.04*(1.+ (10.9491/18.0153)*ovap))
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: LMDZ6/branches/cirrus/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90 (revision 4951)
@@ -49,5 +49,6 @@
prw_ancien, u10m,v10m, treedrg, u_ancien, v_ancien, wake_delta_pbl_TKE, wake_dens, &
ale_bl, ale_bl_trig, alp_bl, &
- ale_wake, ale_bl_stat, AWAKE_S
+ ale_wake, ale_bl_stat, AWAKE_S, &
+ cf_ancien, rvc_ancien
USE comconst_mod, ONLY: pi, dtvr
@@ -94,5 +95,5 @@
USE init_ssrf_m, ONLY: start_init_subsurf
USE phys_state_var_mod, ONLY: beta_aridity, delta_tsurf, awake_dens, cv_gen, &
- ratqs_inter_, rneb_ancien
+ ratqs_inter_
!use ioipsl_getincom
IMPLICIT NONE
@@ -242,4 +243,7 @@
ale_wake=0.
ale_bl_stat=0.
+
+ cf_ancien = 0.
+ rvc_ancien = 0.
z0m(:,:)=0 ! ym missing 5th subsurface initialization
@@ -287,5 +291,4 @@
ratqs_inter_ = 0.002
- rneb_ancien = 0.
CALL phyredem( "startphy.nc" )
Index: LMDZ6/branches/cirrus/libf/misc/readTracFiles_mod.f90
===================================================================
--- LMDZ6/branches/cirrus/libf/misc/readTracFiles_mod.f90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/misc/readTracFiles_mod.f90 (revision 4951)
@@ -114,9 +114,9 @@
!--- SOME PARAMETERS THAT ARE NOT LIKELY TO CHANGE OFTEN
CHARACTER(LEN=maxlen), SAVE :: tran0 = 'air' !--- Default transporting fluid
- CHARACTER(LEN=maxlen), PARAMETER :: old_phases = 'vlirb' !--- Old phases for water (no separator)
- CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsrb' !--- Known phases initials
+ CHARACTER(LEN=maxlen), PARAMETER :: old_phases = 'vlifbc' !--- Old phases for water (no separator)
+ CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsfbc' !--- Known phases initials
INTEGER, PARAMETER :: nphases = LEN_TRIM(known_phases) !--- Number of phases
CHARACTER(LEN=maxlen), SAVE :: phases_names(nphases) & !--- Known phases names
- = ['gaseous', 'liquid ', 'solid ', 'cloud ','blosno ']
+ = ['gaseous', 'liquid ', 'solid ', 'fracld ', 'blosnow', 'cldvapr']
CHARACTER(LEN=1), SAVE :: phases_sep = '_' !--- Phase separator
LOGICAL, SAVE :: tracs_merge = .TRUE. !--- Merge/stack tracers lists
Index: LMDZ6/branches/cirrus/libf/phylmd/clesphys.h
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/clesphys.h (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/clesphys.h (revision 4951)
@@ -94,5 +94,5 @@
LOGICAL :: ok_airs
INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo, NSW, iflag_albedo
- LOGICAL :: ok_ice_sursat, ok_plane_h2o, ok_plane_contrail
+ LOGICAL :: ok_ice_supersat, ok_plane_h2o, ok_plane_contrail
LOGICAL :: ok_chlorophyll
LOGICAL :: ok_strato
@@ -154,5 +154,5 @@
& , ok_lic_melt, ok_lic_cond, aer_type &
& , iflag_rrtm, ok_strato,ok_hines, ok_qch4 &
- & , iflag_ice_thermo, ok_ice_sursat &
+ & , iflag_ice_thermo, ok_ice_supersat &
& , ok_plane_h2o, ok_plane_contrail &
& , ok_gwd_rando, NSW, iflag_albedo &
Index: LMDZ6/branches/cirrus/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/conf_phys_m.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/conf_phys_m.F90 (revision 4951)
@@ -173,5 +173,5 @@
INTEGER,SAVE :: iflag_clw_omp
INTEGER,SAVE :: iflag_ice_thermo_omp
- LOGICAL,SAVE :: ok_ice_sursat_omp
+ LOGICAL,SAVE :: ok_ice_supersat_omp
LOGICAL,SAVE :: ok_plane_h2o_omp, ok_plane_contrail_omp
REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
@@ -1347,15 +1347,15 @@
!
- !Config Key = ok_ice_sursat
- !Config Desc =
- !Config Def = 0
- !Config Help =
- !
- ok_ice_sursat_omp = 0
- CALL getin('ok_ice_sursat',ok_ice_sursat_omp)
+ !Config Key = ok_ice_supersat
+ !Config Desc = include ice supersaturation for cold clouds
+ !Config Def = .FALSE.
+ !Config Help =
+ !
+ ok_ice_supersat_omp = .FALSE.
+ CALL getin('ok_ice_supersat',ok_ice_supersat_omp)
!Config Key = ok_plane_h2o
- !Config Desc =
- !Config Def = 0
+ !Config Desc = include H2O emissions from aviation
+ !Config Def = .FALSE.
!Config Help =
!
@@ -1364,6 +1364,6 @@
!Config Key = ok_plane_contrail
- !Config Desc =
- !Config Def = 0
+ !Config Desc = include the formation of contrail cirrus clouds
+ !Config Def = .FALSE.
!Config Help =
!
@@ -2345,5 +2345,5 @@
rad_chau2 = rad_chau2_omp
iflag_ice_thermo = iflag_ice_thermo_omp
- ok_ice_sursat = ok_ice_sursat_omp
+ ok_ice_supersat = ok_ice_supersat_omp
ok_plane_h2o = ok_plane_h2o_omp
ok_plane_contrail = ok_plane_contrail_omp
@@ -2770,5 +2770,5 @@
WRITE(lunout,*) ' rad_chau2 = ',rad_chau2
WRITE(lunout,*) ' iflag_ice_thermo = ',iflag_ice_thermo
- WRITE(lunout,*) ' ok_ice_sursat = ',ok_ice_sursat
+ WRITE(lunout,*) ' ok_ice_supersat = ',ok_ice_supersat
WRITE(lunout,*) ' ok_plane_h2o = ',ok_plane_h2o
WRITE(lunout,*) ' ok_plane_contrail = ',ok_plane_contrail
Index: LMDZ6/branches/cirrus/libf/phylmd/create_etat0_unstruct_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/create_etat0_unstruct_mod.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/create_etat0_unstruct_mod.F90 (revision 4951)
@@ -258,4 +258,8 @@
prw_ancien = 0.
agesno = 0.
+
+ ! LSCP condensation and ice supersaturation
+ cf_ancien = 0.
+ rvc_ancien = 0.
wake_delta_pbl_TKE(:,:,:)=0
@@ -310,5 +314,4 @@
END IF
ratqs_inter_ = 0.002
- rneb_ancien = 0.
CALL gather_omp(cell_area,cell_area_mpi)
Index: LMDZ6/branches/cirrus/libf/phylmd/ice_sursat_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/ice_sursat_mod.F90 (revision 4950)
+++ (revision )
@@ -1,906 +1,0 @@
-MODULE ice_sursat_mod
-
-IMPLICIT NONE
-
-!--flight inventories
-!
-REAL, SAVE, ALLOCATABLE :: flight_m(:,:) !--flown distance m s-1 per cell
-!$OMP THREADPRIVATE(flight_m)
-REAL, SAVE, ALLOCATABLE :: flight_h2o(:,:) !--emitted kg H2O s-1 per cell
-!$OMP THREADPRIVATE(flight_h2o)
-!
-!--Fixed Parameters
-!
-!--safety parameters for ERF function
-REAL, PARAMETER :: erf_lim = 5., eps = 1.e-10
-!
-!--Tuning parameters (and their default values)
-!
-!--chi gère la répartition statistique de la longueur des frontières
-! entre les zones nuages et ISSR/ciel clair sous-saturé. Gamme de valeur :
-! chi > 1, je n'ai pas regardé de limite max (pour chi = 1, la longueur de
-! la frontière entre ne nuage et l'ISSR est proportionnelle à la
-! répartition ISSR/ciel clair sous-sat dans la maille, i.e. il n'y a pas
-! de favorisation de la localisation de l'ISSR près de nuage. Pour chi = inf,
-! le nuage n'est en contact qu'avec de l'ISSR, quelle que soit la taille
-! de l'ISSR dans la maille.)
-!
-!--l_turb est la longueur de mélange pour la turbulence.
-! dans les tests, ça n'a jamais été modifié pour l'instant.
-!
-!--tun_N est le paramètre qui contrôle l'importance relative de N_2 par rapport à N_1.
-! La valeur est comprise entre 1 et 2 (tun_N = 1 => N_1 = N_2)
-!
-!--tun_ratqs : paramètre qui modifie ratqs en fonction de la valeur de
-! alpha_cld selon la formule ratqs_new = ratqs_old / ( 1 + tun_ratqs *
-! alpha_cld ). Dans le rapport il est appelé beta. Il varie entre 0 et 5
-! (tun_ratqs = 0 => pas de modification de ratqs).
-!
-!--gamma0 and Tgamma: define RHcrit limit above which heterogeneous freezing occurs as a function of T
-!--Karcher and Lohmann (2002)
-!--gamma = 2.583 - t / 207.83
-!--Ren and MacKenzie (2005) reused by Kärcher
-!--gamma = 2.349 - t / 259.0
-!
-!--N_cld: number of clouds in cell (needs to be parametrized at some point)
-!
-!--contrail cross section: typical value found in Freudenthaler et al, GRL, 22, 3501-3504, 1995
-!--in m2, 1000x200 = 200 000 m2 after 15 min
-!
-REAL, SAVE :: chi=1.1, l_turb=50.0, tun_N=1.3, tun_ratqs=3.0
-REAL, SAVE :: gamma0=2.349, Tgamma=259.0, N_cld=100, contrail_cross_section=200000.0
-!$OMP THREADPRIVATE(chi,l_turb,tun_N,tun_ratqs,gamma0,Tgamma,N_cld,contrail_cross_section)
-
-CONTAINS
-
-!*******************************************************************
-!
-SUBROUTINE ice_sursat_init()
-
- USE print_control_mod, ONLY: lunout
- USE ioipsl_getin_p_mod, ONLY : getin_p
-
- IMPLICIT NONE
-
- CALL getin_p('flag_chi',chi)
- CALL getin_p('flag_l_turb',l_turb)
- CALL getin_p('flag_tun_N',tun_N)
- CALL getin_p('flag_tun_ratqs',tun_ratqs)
- CALL getin_p('gamma0',gamma0)
- CALL getin_p('Tgamma',Tgamma)
- CALL getin_p('N_cld',N_cld)
- CALL getin_p('contrail_cross_section',contrail_cross_section)
-
- WRITE(lunout,*) 'Parameters for ice_sursat param'
- WRITE(lunout,*) 'flag_chi = ', chi
- WRITE(lunout,*) 'flag_l_turb = ', l_turb
- WRITE(lunout,*) 'flag_tun_N = ', tun_N
- WRITE(lunout,*) 'flag_tun_ratqs = ', tun_ratqs
- WRITE(lunout,*) 'gamma0 = ', gamma0
- WRITE(lunout,*) 'Tgamma = ', Tgamma
- WRITE(lunout,*) 'N_cld = ', N_cld
- WRITE(lunout,*) 'contrail_cross_section = ', contrail_cross_section
-
-END SUBROUTINE ice_sursat_init
-
-!*******************************************************************
-!
-SUBROUTINE airplane(debut,pphis,pplay,paprs,t_seri)
-
- USE dimphy
- USE mod_grid_phy_lmdz, ONLY: klon_glo
- USE geometry_mod, ONLY: cell_area
- USE phys_cal_mod, ONLY : mth_cur
- USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
- USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
- USE mod_phys_lmdz_para, ONLY: scatter, bcast
- USE print_control_mod, ONLY: lunout
-
- IMPLICIT NONE
-
- INCLUDE "YOMCST.h"
- INCLUDE 'netcdf.inc'
-
- !--------------------------------------------------------
- !--input variables
- !--------------------------------------------------------
- LOGICAL, INTENT(IN) :: debut
- REAL, INTENT(IN) :: pphis(klon), pplay(klon,klev), paprs(klon,klev+1), t_seri(klon,klev)
-
- !--------------------------------------------------------
- ! ... Local variables
- !--------------------------------------------------------
-
- CHARACTER (LEN=20) :: modname='airplane_mod'
- INTEGER :: i, k, kori, iret, varid, error, ncida, klona
- INTEGER,SAVE :: nleva, ntimea
-!$OMP THREADPRIVATE(nleva,ntimea)
- REAL, ALLOCATABLE :: pkm_airpl_glo(:,:,:) !--km/s
- REAL, ALLOCATABLE :: ph2o_airpl_glo(:,:,:) !--molec H2O/cm3/s
- REAL, ALLOCATABLE, SAVE :: zmida(:), zinta(:)
- REAL, ALLOCATABLE, SAVE :: pkm_airpl(:,:,:)
- REAL, ALLOCATABLE, SAVE :: ph2o_airpl(:,:,:)
-!$OMP THREADPRIVATE(pkm_airpl,ph2o_airpl,zmida,zinta)
- REAL :: zalt(klon,klev+1)
- REAL :: zrho, zdz(klon,klev), zfrac
-
- !
- IF (debut) THEN
- !--------------------------------------------------------------------------------
- ! ... Open the file and read airplane emissions
- !--------------------------------------------------------------------------------
- !
- IF (is_mpi_root .AND. is_omp_root) THEN
- !
- iret = nf_open('aircraft_phy.nc', 0, ncida)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to open aircraft_phy.nc file',1)
- ! ... Get lengths
- iret = nf_inq_dimid(ncida, 'time', varid)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimid in aircraft_phy.nc file',1)
- iret = nf_inq_dimlen(ncida, varid, ntimea)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimlen aircraft_phy.nc file',1)
- iret = nf_inq_dimid(ncida, 'vector', varid)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimid aircraft_phy.nc file',1)
- iret = nf_inq_dimlen(ncida, varid, klona)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimlen aircraft_phy.nc file',1)
- iret = nf_inq_dimid(ncida, 'lev', varid)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
- iret = nf_inq_dimlen(ncida, varid, nleva)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimlen aircraft_phy.nc file',1)
- !
- IF ( klona /= klon_glo ) THEN
- WRITE(lunout,*) 'klona & klon_glo =', klona, klon_glo
- CALL abort_physic(modname,'problem klon in aircraft_phy.nc file',1)
- ENDIF
- !
- IF ( ntimea /= 12 ) THEN
- WRITE(lunout,*) 'ntimea=', ntimea
- CALL abort_physic(modname,'problem ntime<>12 in aircraft_phy.nc file',1)
- ENDIF
- !
- ALLOCATE(zmida(nleva), STAT=error)
- IF (error /= 0) CALL abort_physic(modname,'problem to allocate zmida',1)
- ALLOCATE(pkm_airpl_glo(klona,nleva,ntimea), STAT=error)
- IF (error /= 0) CALL abort_physic(modname,'problem to allocate pkm_airpl_glo',1)
- ALLOCATE(ph2o_airpl_glo(klona,nleva,ntimea), STAT=error)
- IF (error /= 0) CALL abort_physic(modname,'problem to allocate ph2o_airpl_glo',1)
- !
- iret = nf_inq_varid(ncida, 'lev', varid)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
- iret = nf_get_var_double(ncida, varid, zmida)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read zmida file',1)
- !
- iret = nf_inq_varid(ncida, 'emi_co2_aircraft', varid) !--CO2 as a proxy for m flown -
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_distance dimid aircraft_phy.nc file',1)
- iret = nf_get_var_double(ncida, varid, pkm_airpl_glo)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read pkm_airpl file',1)
- !
- iret = nf_inq_varid(ncida, 'emi_h2o_aircraft', varid)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_h2o_aircraft dimid aircraft_phy.nc file',1)
- iret = nf_get_var_double(ncida, varid, ph2o_airpl_glo)
- IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read ph2o_airpl file',1)
- !
- ENDIF !--is_mpi_root and is_omp_root
- !
- CALL bcast(nleva)
- CALL bcast(ntimea)
- !
- IF (.NOT.ALLOCATED(zmida)) ALLOCATE(zmida(nleva), STAT=error)
- IF (.NOT.ALLOCATED(zinta)) ALLOCATE(zinta(nleva+1), STAT=error)
- !
- ALLOCATE(pkm_airpl(klon,nleva,ntimea))
- ALLOCATE(ph2o_airpl(klon,nleva,ntimea))
- !
- ALLOCATE(flight_m(klon,klev))
- ALLOCATE(flight_h2o(klon,klev))
- !
- CALL bcast(zmida)
- zinta(1)=0.0 !--surface
- DO k=2, nleva
- zinta(k) = (zmida(k-1)+zmida(k))/2.0*1000.0 !--conversion from km to m
- ENDDO
- zinta(nleva+1)=zinta(nleva)+(zmida(nleva)-zmida(nleva-1))*1000.0 !--extrapolation for last interface
- !print *,'zinta=', zinta
- !
- CALL scatter(pkm_airpl_glo,pkm_airpl)
- CALL scatter(ph2o_airpl_glo,ph2o_airpl)
- !
-!$OMP MASTER
- IF (is_mpi_root .AND. is_omp_root) THEN
- DEALLOCATE(pkm_airpl_glo)
- DEALLOCATE(ph2o_airpl_glo)
- ENDIF !--is_mpi_root
-!$OMP END MASTER
-
- ENDIF !--debut
-!
-!--compute altitude of model level interfaces
-!
- DO i = 1, klon
- zalt(i,1)=pphis(i)/RG !--in m
- ENDDO
-!
- DO k=1, klev
- DO i = 1, klon
- zrho=pplay(i,k)/t_seri(i,k)/RD
- zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho/RG
- zalt(i,k+1)=zalt(i,k)+zdz(i,k) !--in m
- ENDDO
- ENDDO
-!
-!--vertical reprojection
-!
- flight_m(:,:)=0.0
- flight_h2o(:,:)=0.0
-!
- DO k=1, klev
- DO kori=1, nleva
- DO i=1, klon
- !--fraction of layer kori included in layer k
- zfrac=max(0.0,min(zalt(i,k+1),zinta(kori+1))-max(zalt(i,k),zinta(kori)))/(zinta(kori+1)-zinta(kori))
- !--reproject
- flight_m(i,k)=flight_m(i,k) + pkm_airpl(i,kori,mth_cur)*zfrac
- !--reproject
- flight_h2o(i,k)=flight_h2o(i,k) + ph2o_airpl(i,kori,mth_cur)*zfrac
- ENDDO
- ENDDO
- ENDDO
-!
- DO k=1, klev
- DO i=1, klon
- !--molec.cm-3.s-1 / (molec/mol) * kg CO2/mol * m2 * m * cm3/m3 / (kg CO2/m) => m s-1 per cell
- flight_m(i,k)=flight_m(i,k)/RNAVO*44.e-3*cell_area(i)*zdz(i,k)*1.e6/16.37e-3
- flight_m(i,k)=flight_m(i,k)*100.0 !--x100 to augment signal to noise
- !--molec.cm-3.s-1 / (molec/mol) * kg H2O/mol * m2 * m * cm3/m3 => kg H2O s-1 per cell
- flight_h2o(i,k)=flight_h2o(i,k)/RNAVO*18.e-3*cell_area(i)*zdz(i,k)*1.e6
- ENDDO
- ENDDO
-!
-END SUBROUTINE airplane
-
-!********************************************************************
-! simple routine to initialise flight_m and test a flight corridor
-!--Olivier Boucher - 2021
-!
-SUBROUTINE flight_init()
- USE dimphy
- USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
- IMPLICIT NONE
- INTEGER :: i
-
- ALLOCATE(flight_m(klon,klev))
- ALLOCATE(flight_h2o(klon,klev))
- !
- flight_m(:,:) = 0.0 !--initialisation
- flight_h2o(:,:) = 0.0 !--initialisation
- !
- DO i=1, klon
- IF (latitude_deg(i).GE.42.0.AND.latitude_deg(i).LE.48.0) THEN
- flight_m(i,38) = 50000.0 !--5000 m of flight/second in grid cell x 10 scaling
- ENDIF
- ENDDO
-
- RETURN
-END SUBROUTINE flight_init
-
-!*******************************************************************
-!--Routine to deal with ice supersaturation
-!--Determines the respective fractions of unsaturated clear sky, ice supersaturated clear sky and cloudy sky
-!--Diagnoses regions prone for non-persistent and persistent contrail formation
-!
-!--Audran Borella - 2021
-!
-SUBROUTINE ice_sursat(pplay, dpaprs, dtime, i, k, t, q, gamma_ss, &
- qsat, t_actuel, rneb_seri, ratqs, rneb, qincld, &
- rnebss, qss, Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
- !
- USE dimphy
- USE print_control_mod, ONLY: prt_level, lunout
- USE phys_state_var_mod, ONLY: pbl_tke, t_ancien
- USE phys_local_var_mod, ONLY: N1_ss, N2_ss
- USE phys_local_var_mod, ONLY: drneb_sub, drneb_con, drneb_tur, drneb_avi
-!! USE phys_local_var_mod, ONLY: Tcontr, qcontr, fcontrN, fcontrP
- USE indice_sol_mod, ONLY: is_ave
- USE geometry_mod, ONLY: cell_area
- !
- IMPLICIT NONE
- INCLUDE "YOMCST.h"
- INCLUDE "YOETHF.h"
- INCLUDE "FCTTRE.h"
- INCLUDE "clesphys.h"
-
- !
- ! Input
- ! Beware: this routine works on a gridpoint!
- !
- REAL, INTENT(IN) :: pplay ! layer pressure (Pa)
- REAL, INTENT(IN) :: dpaprs ! layer delta pressure (Pa)
- REAL, INTENT(IN) :: dtime ! intervalle du temps (s)
- REAL, INTENT(IN) :: t ! température advectée (K)
- REAL, INTENT(IN) :: qsat ! vapeur de saturation
- REAL, INTENT(IN) :: t_actuel ! temperature actuelle de la maille (K)
- REAL, INTENT(IN) :: rneb_seri ! fraction nuageuse en memoire
- INTEGER, INTENT(IN) :: i, k
- !
- ! Input/output
- !
- REAL, INTENT(INOUT) :: q ! vapeur de la maille (=zq)
- REAL, INTENT(INOUT) :: ratqs ! determine la largeur de distribution de vapeur
- REAL, INTENT(INOUT) :: Tcontr, qcontr, qcontr2, fcontrN, fcontrP
- !
- ! Output
- !
- REAL, INTENT(OUT) :: gamma_ss !
- REAL, INTENT(OUT) :: rneb ! cloud fraction
- REAL, INTENT(OUT) :: qincld ! in-cloud total water
- REAL, INTENT(OUT) :: rnebss ! ISSR fraction
- REAL, INTENT(OUT) :: qss ! in-ISSR total water
- !
- ! Local
- !
- REAL PI
- PARAMETER (PI=4.*ATAN(1.))
- REAL rnebclr, gamma_prec
- REAL qclr, qvc, qcld, qi
- REAL zrho, zdz, zrhodz
- REAL pdf_N, pdf_N1, pdf_N2
- REAL pdf_a, pdf_b
- REAL pdf_e1, pdf_e2, pdf_k
- REAL drnebss, drnebclr, dqss, dqclr, sum_rneb_rnebss, dqss_avi
- REAL V_cell !--volume of the cell
- REAL M_cell !--dry mass of the cell
- REAL tke, sig, L_tur, b_tur, q_eq
- REAL V_env, V_cld, V_ss, V_clr
- REAL zcor
- !
- !--more local variables for diagnostics
- !--imported from YOMCST.h
- !--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1)
- !--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air
- !--values from Schumann, Meteorol Zeitschrift, 1996
- !--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
- !--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
- REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1)
- REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1)
- REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft
- !--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
- !--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
- REAL :: Gcontr
- !--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2)
- !--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1)
- REAL :: qsatliqcontr
-
- ! Initialisations
- zrho = pplay / t / RD !--dry density kg m-3
- zrhodz = dpaprs / RG !--dry air mass kg m-2
- zdz = zrhodz / zrho !--cell thickness m
- V_cell = zdz * cell_area(i) !--cell volume m3
- M_cell = zrhodz * cell_area(i) !--cell dry air mass kg
- !
- ! Recuperation de la memoire sur la couverture nuageuse
- rneb = rneb_seri
- !
- ! Ajout des émissions de H2O dues à l'aviation
- ! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
- ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O )
- ! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
- ! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
- ! flight_h2O is in kg H2O / s / cell
- !
- IF (ok_plane_h2o) THEN
- q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) )
- ENDIF
- !
- !--Estimating gamma
- gamma_ss = MAX(1.0, gamma0 - t_actuel/Tgamma)
- !gamma_prec = MAX(1.0, gamma0 - t_ancien(i,k)/Tgamma) !--formulation initiale d Audran
- gamma_prec = MAX(1.0, gamma0 - t/Tgamma) !--autre formulation possible basée sur le T du pas de temps
- !
- ! Initialisation de qvc : q_sat du pas de temps precedent
- !qvc = R2ES*FOEEW(t_ancien(i,k),1.)/pplay !--formulation initiale d Audran
- qvc = R2ES*FOEEW(t,1.)/pplay !--autre formulation possible basée sur le T du pas de temps
- qvc = min(0.5,qvc)
- zcor = 1./(1.-RETV*qvc)
- qvc = qvc*zcor
- !
- ! Modification de ratqs selon formule proposee : ksi_new = ksi_old/(1+beta*alpha_cld)
- ratqs = ratqs / (tun_ratqs*rneb_seri + 1.)
- !
- ! Calcul de N
- pdf_k = -sqrt(log(1.+ratqs**2.))
- pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
- pdf_b = pdf_k/(2.*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- pdf_N = max(0.,sign(rneb,pdf_e1))
- ELSE
- pdf_e1 = erf(pdf_e1)
- pdf_e1 = 0.5*(1.+pdf_e1)
- pdf_N = max(0.,rneb/pdf_e1)
- ENDIF
- !
- ! On calcule ensuite N1 et N2. Il y a deux cas : N > 1 et N <= 1
- ! Cas 1 : N > 1. N'arrive en theorie jamais, c'est une barriere
- ! On perd la memoire sur la temperature (sur qvc) pour garder
- ! celle sur alpha_cld
- IF (pdf_N.GT.1.) THEN
- ! On inverse alpha_cld = int_qvc^infty P(q) dq
- ! pour determiner qvc = f(alpha_cld)
- ! On approxime en serie entiere erf-1(x)
- qvc = 2.*rneb-1.
- qvc = qvc + PI/12.*qvc**3 + 7.*PI**2/480.*qvc**5 &
- + 127.*PI**3/40320.*qvc**7 + 4369.*PI**4/5806080.*qvc**9 &
- + 34807.*PI**5/182476800.*qvc**11
- qvc = sqrt(PI)/2.*qvc
- qvc = (qvc-pdf_b)*pdf_k*sqrt(2.)
- qvc = q*exp(qvc)
-
- ! On met a jour rneb avec la nouvelle valeur de qvc
- ! La maj est necessaire a cause de la serie entiere approximative
- pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
- pdf_e1 = 0.5*(1.+pdf_e1)
- rneb = pdf_e1
-
- ! Si N > 1, N1 et N2 = 1
- pdf_N1 = 1.
- pdf_N2 = 1.
-
- ! Cas 2 : N <= 1
- ELSE
- ! D'abord on calcule N2 avec le tuning
- pdf_N2 = min(1.,pdf_N*tun_N)
-
- ! Puis N1 pour assurer la conservation de alpha_cld
- pdf_a = log(qvc*gamma_prec/q)/(pdf_k*sqrt(2.))
- pdf_e2 = pdf_a+pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
- pdf_e2 = 0.5*(1.+pdf_e2) ! integrale sous P pour q > gamma qsat
-
- IF (abs(pdf_e1-pdf_e2).LT.eps) THEN
- pdf_N1 = pdf_N2
- ELSE
- pdf_N1 = (rneb-pdf_N2*pdf_e2)/(pdf_e1-pdf_e2)
- ENDIF
-
- ! Barriere qui traite le cas gamma_prec = 1.
- IF (pdf_N1.LE.0.) THEN
- pdf_N1 = 0.
- IF (pdf_e2.GT.eps) THEN
- pdf_N2 = rneb/pdf_e2
- ELSE
- pdf_N2 = 0.
- ENDIF
- ENDIF
- ENDIF
-
- ! Physique 1
- ! Sublimation
- IF (qvc.LT.qsat) THEN
- pdf_a = log(qvc/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
-
- pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
- pdf_e2 = pdf_a+pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
-
- pdf_e1 = 0.5*pdf_N1*(pdf_e1-pdf_e2)
-
- ! Calcul et ajout de la tendance
- drneb_sub(i,k) = - pdf_e1/dtime !--unit [s-1]
- rneb = rneb + drneb_sub(i,k)*dtime
- ELSE
- drneb_sub(i,k) = 0.
- ENDIF
-
- ! NOTE : verifier si ca marche bien pour gamma proche de 1.
-
- ! Condensation
- IF (gamma_ss*qsat.LT.gamma_prec*qvc) THEN
-
- pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
-
- pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.))
- pdf_e2 = pdf_a+pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
-
- pdf_e1 = 0.5*(1.-pdf_N1)*(pdf_e1-pdf_e2)
- pdf_e2 = 0.5*(1.-pdf_N2)*(1.+pdf_e2)
-
- ! Calcul et ajout de la tendance
- drneb_con(i,k) = (pdf_e1 + pdf_e2)/dtime !--unit [s-1]
- rneb = rneb + drneb_con(i,k)*dtime
-
- ELSE
-
- pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
- pdf_e1 = 0.5*(1.-pdf_N2)*(1.+pdf_e1)
-
- ! Calcul et ajout de la tendance
- drneb_con(i,k) = pdf_e1/dtime !--unit [s-1]
- rneb = rneb + drneb_con(i,k)*dtime
-
- ENDIF
-
- ! Calcul des grandeurs diagnostiques
- ! Determination des grandeurs ciel clair
- pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
- pdf_e1 = 0.5*(1.-pdf_e1)
-
- pdf_e2 = pdf_a-pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
- pdf_e2 = 0.5*q*(1.-pdf_e2)
-
- rnebclr = pdf_e1
- qclr = pdf_e2
-
- ! Determination de q_cld
- ! Partie 1
- pdf_a = log(max(qsat,qvc)/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a-pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
-
- pdf_a = log(min(gamma_ss*qsat,gamma_prec*qvc)/q)/(pdf_k*sqrt(2.))
- pdf_e2 = pdf_a-pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
-
- pdf_e1 = 0.5*q*pdf_N1*(pdf_e1-pdf_e2)
-
- qcld = pdf_e1
-
- ! Partie 2 (sous condition)
- IF (gamma_ss*qsat.GT.gamma_prec*qvc) THEN
- pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a-pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
-
- pdf_e2 = 0.5*q*pdf_N2*(pdf_e2-pdf_e1)
-
- qcld = qcld + pdf_e2
-
- pdf_e2 = pdf_e1
- ENDIF
-
- ! Partie 3
- pdf_e2 = 0.5*q*(1.+pdf_e2)
-
- qcld = qcld + pdf_e2
-
- ! Fin du calcul de q_cld
-
- ! Determination des grandeurs ISSR via les equations de conservation
- rneb=MIN(rneb, 1. - rnebclr - eps) !--ajout OB - barrière
- rnebss = MAX(0.0, 1. - rnebclr - rneb) !--ajout OB
- qss = MAX(0.0, q - qclr - qcld) !--ajout OB
-
- ! Physique 2 : Turbulence
- IF (rneb.GT.eps.AND.rneb.LT.1.-eps) THEN ! rneb != 0 and != 1
- !
- tke = pbl_tke(i,k,is_ave)
- ! A MODIFIER formule a revoir
- L_tur = min(l_turb, sqrt(tke)*dtime)
-
- ! On fait pour l'instant l'hypothese a = 3b. V = 4/3 pi a b**2 = alpha_cld
- ! donc b = alpha_cld/4pi **1/3.
- b_tur = (rneb*V_cell/4./PI/N_cld)**(1./3.)
- ! On verifie que la longeur de melange n'est pas trop grande
- IF (L_tur.GT.b_tur) THEN
- L_tur = b_tur
- ENDIF
-
- V_env = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. + 3.*b_tur*(L_tur**2.))
- V_cld = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. - 3.*b_tur*(L_tur**2.))
- V_cld = abs(V_cld)
-
- ! Repartition statistique des zones de contact nuage-ISSR et nuage-ciel clair
- sig = rnebss/(chi*rnebclr+rnebss)
- V_ss = MIN(sig*V_env,rnebss*V_cell)
- V_clr = MIN((1.-sig)*V_env,rnebclr*V_cell)
- V_cld = MIN(V_cld,rneb*V_cell)
- V_env = V_ss + V_clr
-
- ! ISSR => rneb
- drnebss = -1.*V_ss/V_cell
- dqss = drnebss*qss/MAX(eps,rnebss)
-
- ! Clear sky <=> rneb
- q_eq = V_env*qclr/MAX(eps,rnebclr) + V_cld*qcld/MAX(eps,rneb)
- q_eq = q_eq/(V_env + V_cld)
-
- IF (q_eq.GT.qsat) THEN
- drnebclr = - V_clr/V_cell
- dqclr = drnebclr*qclr/MAX(eps,rnebclr)
- ELSE
- drnebclr = V_cld/V_cell
- dqclr = drnebclr*qcld/MAX(eps,rneb)
- ENDIF
-
- ! Maj des variables avec les tendances
- rnebclr = MAX(0.0,rnebclr + drnebclr) !--OB ajout d'un max avec eps (il faudrait modified drnebclr pour le diag)
- qclr = MAX(0.0, qclr + dqclr) !--OB ajout d'un max avec 0
-
- rneb = rneb - drnebclr - drnebss
- qcld = qcld - dqclr - dqss
-
- rnebss = MAX(0.0,rnebss + drnebss) !--OB ajout d'un max avec eps (il faudrait modifier drnebss pour le diag)
- qss = MAX(0.0, qss + dqss) !--OB ajout d'un max avec 0
-
- ! Tendances pour le diagnostic
- drneb_tur(i,k) = (drnebclr + drnebss)/dtime !--unit [s-1]
-
- ENDIF ! rneb
-
- !--add a source of cirrus from aviation contrails
- IF (ok_plane_contrail) THEN
- drneb_avi(i,k) = rnebss*flight_m(i,k)*contrail_cross_section/V_cell !--tendency rneb due to aviation [s-1]
- drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime) !--majoration
- dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss) !--tendency q aviation [kg kg-1 s-1]
- rneb = rneb + drneb_avi(i,k)*dtime !--add tendency to rneb
- qcld = qcld + dqss_avi*dtime !--add tendency to qcld
- rnebss = rnebss - drneb_avi(i,k)*dtime !--add tendency to rnebss
- qss = qss - dqss_avi*dtime !--add tendency to qss
- ELSE
- drneb_avi(i,k)=0.0
- ENDIF
-
- ! Barrieres
- ! ISSR trop petite
- IF (rnebss.LT.eps) THEN
- rneb = MIN(rneb + rnebss,1.0-eps) !--ajout OB barriere
- qcld = qcld + qss
- rnebss = 0.
- qss = 0.
- ENDIF
-
- ! le nuage est trop petit
- IF (rneb.LT.eps) THEN
- ! s'il y a une ISSR on met tout dans l'ISSR, sinon dans le
- ! clear sky
- IF (rnebss.LT.eps) THEN
- rnebclr = 1.
- rnebss = 0. !--ajout OB
- qclr = q
- ELSE
- rnebss = MIN(rnebss + rneb,1.0-eps) !--ajout OB barriere
- qss = qss + qcld
- ENDIF
- rneb = 0.
- qcld = 0.
- qincld = qsat * gamma_ss
- ELSE
- qincld = qcld / rneb
- ENDIF
-
- !--OB ajout borne superieure
- sum_rneb_rnebss=rneb+rnebss
- rneb=rneb*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss)
- rnebss=rnebss*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss)
-
- ! On ecrit dans la memoire
- N1_ss(i,k) = pdf_N1
- N2_ss(i,k) = pdf_N2
-
- !--Diagnostics only used from last iteration
- !--test
- !!Tcontr(i,k)=200.
- !!fcontrN(i,k)=1.0
- !!fcontrP(i,k)=0.5
- !
- !--slope of dilution line in exhaust
- !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
- Gcontr = EiH2O * RCPD * pplay / (eps_w*Qheat*(1.-eta)) !--Pa K-1
- !--critical T_LM below which no liquid contrail can form in exhaust
- !Tcontr(i,k) = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K
- IF (Gcontr .GT. 0.1) THEN
- !
- Tcontr = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K
- !print *,'Tcontr=',iter,i,k,eps_w,pplay,Gcontr,Tcontr(i,k)
- !--Psat with index 0 in FOEEW to get saturation wrt liquid water corresponding to Tcontr
- !qsatliqcontr = RESTT*FOEEW(Tcontr(i,k),0.) !--Pa
- qsatliqcontr = RESTT*FOEEW(Tcontr,0.) !--Pa
- !--Critical water vapour above which there is contrail formation for ambiant temperature
- !qcontr(i,k) = Gcontr*(t-Tcontr(i,k)) + qsatliqcontr !--Pa
- qcontr = Gcontr*(t-Tcontr) + qsatliqcontr !--Pa
- !--Conversion of qcontr in specific humidity - method 1
- !qcontr(i,k) = RD/RV*qcontr(i,k)/pplay !--so as to return to something similar to R2ES*FOEEW/pplay
- qcontr2 = RD/RV*qcontr/pplay !--so as to return to something similar to R2ES*FOEEW/pplay
- !qcontr(i,k) = min(0.5,qcontr(i,k)) !--and then we apply the same correction term as for qsat
- qcontr2 = min(0.5,qcontr2) !--and then we apply the same correction term as for qsat
- !zcor = 1./(1.-RETV*qcontr(i,k)) !--for consistency with qsat but is it correct at all?
- zcor = 1./(1.-RETV*qcontr2) !--for consistency with qsat but is it correct at all as p is dry?
- !zcor = 1./(1.+qcontr2) !--for consistency with qsat
- !qcontr(i,k) = qcontr(i,k)*zcor
- qcontr2 = qcontr2*zcor
- qcontr2=MAX(1.e-10,qcontr2) !--eliminate negative values due to extrapolation on dilution curve
- !--Conversion of qcontr in specific humidity - method 2
- !qcontr(i,k) = eps_w*qcontr(i,k) / (pplay+eps_w*qcontr(i,k))
- !qcontr=MAX(1.E-10,qcontr)
- !qcontr2 = eps_w*qcontr / (pplay+eps_w*qcontr)
- !
- IF (t .LT. Tcontr) THEN !--contrail formation is possible
- !
- !--compute fractions of persistent (P) and non-persistent(N) contrail potential regions
- !!IF (qcontr(i,k).GE.qsat) THEN
- IF (qcontr2.GE.qsat) THEN
- !--none of the unsaturated clear sky is prone for contrail formation
- !!fcontrN(i,k) = 0.0
- fcontrN = 0.0
- !
- !--integral of P(q) from qsat to qcontr in ISSR
- pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
- !
- !!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.))
- pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.))
- pdf_e2 = pdf_a+pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
- !
- !!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2))
- fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2))
- !
- pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
- !
- !!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.))
- pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.))
- pdf_e2 = pdf_a+pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
- !
- !!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2))
- fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2))
- !
- pdf_a = log(MAX(qsat,qvc)/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
- !
- !!pdf_a = log(MIN(qcontr(i,k),MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.))
- pdf_a = log(MIN(qcontr2,MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.))
- pdf_e2 = pdf_a+pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
- !
- !!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2))
- fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2))
- !
- pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
- !
- !!pdf_a = log(MIN(qcontr(i,k),gamma_ss*qsat)/q)/(pdf_k*sqrt(2.))
- pdf_a = log(MIN(qcontr2,gamma_ss*qsat)/q)/(pdf_k*sqrt(2.))
- pdf_e2 = pdf_a+pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
- !
- !!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2))
- fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2))
- !
- ELSE !--qcontr LT qsat
- !
- !--all of ISSR is prone for contrail formation
- !!fcontrP(i,k) = rnebss
- fcontrP = rnebss
- !
- !--integral of zq from qcontr to qsat in unsaturated clear-sky region
- !!pdf_a = log(qcontr(i,k)/q)/(pdf_k*sqrt(2.))
- pdf_a = log(qcontr2/q)/(pdf_k*sqrt(2.))
- pdf_e1 = pdf_a+pdf_b !--normalement pdf_b est deja defini
- IF (abs(pdf_e1).GE.erf_lim) THEN
- pdf_e1 = sign(1.,pdf_e1)
- ELSE
- pdf_e1 = erf(pdf_e1)
- ENDIF
- !
- pdf_a = log(qsat/q)/(pdf_k*sqrt(2.))
- pdf_e2 = pdf_a+pdf_b
- IF (abs(pdf_e2).GE.erf_lim) THEN
- pdf_e2 = sign(1.,pdf_e2)
- ELSE
- pdf_e2 = erf(pdf_e2)
- ENDIF
- !
- !!fcontrN(i,k) = 0.5*(pdf_e1-pdf_e2)
- fcontrN = 0.5*(pdf_e1-pdf_e2)
- !!fcontrN=2.0
- !
- ENDIF
- !
- ENDIF !-- t < Tcontr
- !
- ENDIF !-- Gcontr > 0.1
-
- RETURN
-END SUBROUTINE ice_sursat
-!
-!*******************************************************************
-!
-END MODULE ice_sursat_mod
Index: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90 (revision 4951)
@@ -8,5 +8,5 @@
SUBROUTINE lscp(klon,klev,dtime,missing_val, &
paprs,pplay,temp,qt,ptconv,ratqs, &
- d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri, &
+ d_t, d_q, d_ql, d_qi, rneb, rneblsvol, &
pfraclr,pfracld, &
radocond, radicefrac, rain, snow, &
@@ -14,8 +14,12 @@
prfl, psfl, rhcl, qta, fraca, &
tv, pspsk, tla, thl, iflag_cld_th, &
- iflag_ice_thermo, ok_ice_sursat, qsatl, qsats, &
- distcltop,temp_cltop, &
- qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
- Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
+ iflag_ice_thermo, distcltop, temp_cltop, cell_area,&
+ cf_seri, rvc_seri, u_seri, v_seri, pbl_eps, &
+ qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
+ ratio_qi_qtot, dcf_sub, dcf_con, dcf_mix, &
+ dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, &
+ dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati, &
+ Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi,&
+ dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, &
@@ -98,9 +102,9 @@
USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat
USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
-USE ice_sursat_mod, ONLY : ice_sursat
+USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
USE lmdz_lscp_poprecip, ONLY : poprecip_precld, poprecip_postcld
! Use du module lmdz_lscp_ini contenant les constantes
-USE lmdz_lscp_ini, ONLY : prt_level, lunout
+USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
@@ -111,4 +115,5 @@
USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
USE lmdz_lscp_ini, ONLY : ok_poprecip
+USE lmdz_lscp_ini, ONLY : ok_external_lognormal, ok_ice_supersat, ok_unadjusted_clouds
IMPLICIT NONE
@@ -132,6 +137,4 @@
INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics
! CR: if iflag_ice_thermo=2, only convection is active
- LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated
-
LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active
@@ -150,7 +153,18 @@
REAL, DIMENSION(klon,klev), INTENT(INOUT) :: ratqs ! function of pressure that sets the large-scale
-
- ! Input sursaturation en glace
- REAL, DIMENSION(klon,klev), INTENT(INOUT):: rneb_seri ! fraction nuageuse en memoire
+ ! INPUT/OUTPUT condensation and ice supersaturation
+ !--------------------------------------------------
+ REAL, DIMENSION(klon,klev), INTENT(INOUT):: cf_seri ! cloud fraction [-]
+ REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratio_qi_qtot ! solid specific water to total specific water ratio [-]
+ REAL, DIMENSION(klon,klev), INTENT(INOUT):: rvc_seri ! cloudy water vapor to total water vapor ratio [-]
+ REAL, DIMENSION(klon,klev), INTENT(IN) :: u_seri ! eastward wind [m/s]
+ REAL, DIMENSION(klon,klev), INTENT(IN) :: v_seri ! northward wind [m/s]
+ REAL, DIMENSION(klon,klev+1), INTENT(IN) :: pbl_eps ! TKE dissipation [?]
+ REAL, DIMENSION(klon), INTENT(IN) :: cell_area ! area of each cell [m2]
+
+ ! INPUT/OUTPUT aviation
+ !--------------------------------------------------
+ REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_dist ! Aviation distance flown within the mesh [m/s/mesh]
+ REAL, DIMENSION(klon,klev), INTENT(IN) :: flight_h2o ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
! OUTPUT variables
@@ -170,6 +184,4 @@
REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2]
REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl ! saturation specific humidity wrt liquid [kg/kg]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsats ! saturation specific humidity wrt ice [kg/kg]
REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2]
REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2]
@@ -183,22 +195,37 @@
REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-]
- ! for supersaturation and contrails parameterisation
-
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: qclr ! specific total water content in clear sky region [kg/kg]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld ! specific total water content in cloudy region [kg/kg]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: qss ! specific total water content in supersat region [kg/kg]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: qvc ! specific vapor content in clouds [kg/kg]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebclr ! mesh fraction of clear sky [-]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebss ! mesh fraction of ISSR [-]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcontr ! threshold temperature for contrail formation [K]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr ! threshold humidity for contrail formation [kg/kg]
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q)
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrN ! fraction of grid favourable to non-persistent contrails
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrP ! fraction of grid favourable to persistent contrails
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals
- REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment
+ ! for condensation and ice supersaturation
+
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsub !--specific total water content in sub-saturated clear sky region [kg/kg]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: qissr !--specific total water content in supersat region [kg/kg]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld !--specific total water content in cloudy region [kg/kg]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: subfra !--mesh fraction of subsaturated clear sky [-]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: issrfra !--mesh fraction of ISSR [-]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_cond !--coefficient governing the ice nucleation RHi threshold [-]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_sub !--cloud fraction tendency because of sublimation [s-1]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_con !--cloud fraction tendency because of condensation [s-1]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_mix !--cloud fraction tendency because of cloud mixing [s-1]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_adj !--specific ice content tendency because of temperature adjustment [kg/kg/s]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_sub !--specific ice content tendency because of sublimation [kg/kg/s]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_con !--specific ice content tendency because of condensation [kg/kg/s]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_mix !--specific ice content tendency because of cloud mixing [kg/kg/s]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_adj !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_sub !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_con !--specific cloud water vapor tendency because of condensation [kg/kg/s]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_mix !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl !--saturation specific humidity wrt liquid [kg/kg]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsati !--saturation specific humidity wrt ice [kg/kg]
+
+ ! for contrails and aviation
+
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcontr !--threshold temperature for contrail formation [K]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr !--threshold humidity for contrail formation [kg/kg]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr2 !--// (2nd expression more consistent with LMDZ expression of q)
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrN !--fraction of grid favourable to non-persistent contrails
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrP !--fraction of grid favourable to persistent contrails
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_avi !--cloud fraction tendency because of aviation [s-1]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_avi !--specific ice content tendency because of aviation [kg/kg/s]
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_avi !--specific cloud water vapor tendency because of aviation [kg/kg/s]
+
! for POPRECIP
@@ -217,4 +244,11 @@
REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsfreez !--snow tendency due to freezing [kg/kg/s]
REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrfreez !--rain tendency due to freezing [kg/kg/s]
+
+ ! for thermals
+
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth !--mean saturation deficit in thermals
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv !--mean saturation deficit in environment
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath !--std of saturation deficit in thermals
+ REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv !--std of saturation deficit in environment
@@ -282,4 +316,15 @@
REAL :: effective_zneb
REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D
+
+ ! for condensation and ice supersaturation
+ REAL, DIMENSION(klon) :: qvc, shear
+ REAL :: delta_z
+ !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
+ ! Constants used for calculating ratios that are advected (using a parent-child
+ ! formalism). This is not done in the dynamical core because at this moment,
+ ! only isotopes can use this parent-child formalism. Note that the two constants
+ ! are the same as the one use in the dynamical core, being also defined in
+ ! dyn3d_common/infotrac.F90
+ REAL :: min_qParent, min_ratio
@@ -363,15 +408,36 @@
temp_cltop1D(:) = 0.0
ztupnew(:)=0.0
-!--ice supersaturation
-gamma_ss(:,:) = 1.0
-qss(:,:) = 0.0
-rnebss(:,:) = 0.0
-Tcontr(:,:) = missing_val
-qcontr(:,:) = missing_val
-qcontr2(:,:) = missing_val
-fcontrN(:,:) = 0.0
-fcontrP(:,:) = 0.0
+
distcltop(:,:)=0.
temp_cltop(:,:)=0.
+
+!--Ice supersaturation
+gamma_cond(:,:) = 1.
+qissr(:,:) = 0.
+issrfra(:,:) = 0.
+dcf_sub(:,:) = 0.
+dcf_con(:,:) = 0.
+dcf_mix(:,:) = 0.
+dqi_adj(:,:) = 0.
+dqi_sub(:,:) = 0.
+dqi_con(:,:) = 0.
+dqi_mix(:,:) = 0.
+dqvc_adj(:,:) = 0.
+dqvc_sub(:,:) = 0.
+dqvc_con(:,:) = 0.
+dqvc_mix(:,:) = 0.
+fcontrN(:,:) = 0.
+fcontrP(:,:) = 0.
+Tcontr(:,:) = missing_val
+qcontr(:,:) = missing_val
+qcontr2(:,:) = missing_val
+dcf_avi(:,:) = 0.
+dqi_avi(:,:) = 0.
+dqvc_avi(:,:) = 0.
+qvc(:) = 0.
+shear(:) = 0.
+min_qParent = 1.e-30
+min_ratio = 1.e-16
+
!-- poprecip
qraindiag(:,:)= 0.
@@ -759,4 +825,6 @@
rneblsvol(i,k)=ctot_vol(i,k)
zqn(i)=qcloud(i)
+ !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+ qvc(i) = rneb(i,k) * zqs(i)
ENDDO
@@ -794,12 +862,15 @@
DO iter=1,iflag_fisrtilp_qsat+1
+ keepgoing(:) = .FALSE.
+
DO i=1,klon
! keepgoing = .true. while convergence is not satisfied
- keepgoing(i)=ABS(DT(i)).GT.DDT0
-
- IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
+
+ IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
! if not convergence:
+ ! we calculate a new iteration
+ keepgoing(i) = .TRUE.
! P2.2.1> cloud fraction and condensed water mass calculation
@@ -833,7 +904,65 @@
CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:))
+
+ !--AB Activates a condensation scheme that allows for
+ !--ice supersaturation and contrails evolution from aviation
+ IF (ok_ice_supersat) THEN
+
+ !--Calculate the shear value (input for condensation and ice supersat)
+ DO i = 1, klon
+ !--Cell thickness [m]
+ delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
+ IF ( iftop ) THEN
+ ! top
+ shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
+ + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
+ ELSEIF ( k .EQ. 1 ) THEN
+ ! surface
+ shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
+ + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
+ ELSE
+ ! other layers
+ shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
+ - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
+ + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
+ - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
+ ENDIF
+ ENDDO
+
+ !---------------------------------------------
+ !-- CONDENSATION AND ICE SUPERSATURATION --
+ !---------------------------------------------
+
+ CALL condensation_ice_supersat( &
+ klon, dtime, missing_val, &
+ pplay(:,k), paprs(:,k), paprs(:,k+1), &
+ cf_seri(:,k), rvc_seri(:,k), ratio_qi_qtot(:,k), &
+ shear(:), pbl_eps(:,k), cell_area(:), &
+ Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), keepgoing(:), &
+ rneb(:,k), zqn(:), qvc(:), issrfra(:,k), qissr(:,k), &
+ dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
+ dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
+ dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
+ Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
+ flight_dist(:,k), flight_h2o(:,k), &
+ dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
+
+
+ !--If .TRUE., calls an externalised version of the generalised
+ !--lognormal condensation scheme (Bony and Emanuel 2001)
+ !--Numerically, convergence is conserved with this option
+ !--The objective is to simplify LSCP
+ ELSEIF ( ok_external_lognormal ) THEN
+
+ CALL condensation_lognormal( &
+ klon, Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), &
+ keepgoing(:), rneb(:,k), zqn(:), qvc(:))
+
+
+ ELSE !--Generalised lognormal (Bony and Emanuel 2001)
+
DO i=1,klon !todoan : check if loop in i is needed
- IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
+ IF (keepgoing(i)) THEN
zpdf_sig(i)=ratqs(i,k)*zq(i)
@@ -849,31 +978,23 @@
zpdf_e2(i)=1.-erf(zpdf_e2(i))
- IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN
-
IF (zpdf_e1(i).LT.1.e-10) THEN
rneb(i,k)=0.
- zqn(i)=gammasat(i)*zqs(i)
+ zqn(i)=zqs(i)
+ !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+ qvc(i) = 0.
ELSE
rneb(i,k)=0.5*zpdf_e1(i)
zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
+ !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+ qvc(i) = rneb(i,k) * zqs(i)
ENDIF
- rnebss(i,k)=0.0 !--added by OB (needed because of the convergence loop on time)
- fcontrN(i,k)=0.0 !--idem
- fcontrP(i,k)=0.0 !--idem
- qss(i,k)=0.0 !--idem
-
- ELSE ! in case of ice supersaturation by Audran
-
- !------------------------------------
- ! ICE SUPERSATURATION
- !------------------------------------
-
- CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), &
- gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), &
- rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), &
- Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
-
- ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))
+ ENDIF ! keepgoing
+ ENDDO ! loop on klon
+
+ ENDIF ! .NOT. ok_ice_supersat
+
+ DO i=1,klon
+ IF (keepgoing(i)) THEN
! If vertical heterogeneity, change fraction by volume as well
@@ -958,8 +1079,20 @@
zqn(i) = zq(i)
rneb(i,k) = 1.0
- zcond(i) = MAX(0.0,zqn(i)-zqs(i))
+ IF ( ok_unadjusted_clouds ) THEN
+ !--AB We relax the saturation adjustment assumption
+ !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
+ zcond(i) = MAX(0., zqn(i) - qvc(i))
+ ELSE
+ zcond(i) = MAX(0.0,zqn(i)-zqs(i))
+ ENDIF
rhcl(i,k)=1.0
ELSE
- zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+ IF ( ok_unadjusted_clouds ) THEN
+ !--AB We relax the saturation adjustment assumption
+ !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
+ zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
+ ELSE
+ zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+ ENDIF
! following line is very strange and probably wrong:
rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
@@ -985,4 +1118,47 @@
rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
ENDIF
+
+ !--AB Write diagnostics and tracers for ice supersaturation
+ IF ( ok_ice_supersat ) THEN
+ CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
+ CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
+
+ DO i = 1, klon
+
+ cf_seri(i,k) = rneb(i,k)
+
+ IF ( .NOT. ok_unadjusted_clouds ) THEN
+ qvc(i) = zqs(i) * rneb(i,k)
+ ENDIF
+ IF ( zq(i) .GT. min_qParent ) THEN
+ rvc_seri(i,k) = qvc(i) / zq(i)
+ ELSE
+ rvc_seri(i,k) = min_ratio
+ ENDIF
+ !--The MIN barrier is NEEDED because of:
+ !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
+ !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
+ !--The MAX barrier is a safeguard that should not be activated
+ rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
+
+ !--Diagnostics
+ gamma_cond(i,k) = gammasat(i)
+ IF ( issrfra(i,k) .LT. eps ) THEN
+ issrfra(i,k) = 0.
+ qissr(i,k) = 0.
+ ENDIF
+ subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
+ qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
+ IF ( subfra(i,k) .LT. eps ) THEN
+ subfra(i,k) = 0.
+ qsub(i,k) = 0.
+ ENDIF
+ qcld(i,k) = qvc(i) + zcond(i)
+ IF ( cf_seri(i,k) .LT. eps ) THEN
+ qcld(i,k) = 0.
+ ENDIF
+ ENDDO
+ ENDIF
+
! ----------------------------------------------------------------
@@ -1407,23 +1583,4 @@
ENDDO
- !--save some variables for ice supersaturation
- !
- DO i = 1, klon
- ! for memory
- rneb_seri(i,k) = rneb(i,k)
-
- ! for diagnostics
- rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
-
- qvc(i,k) = zqs(i) * rneb(i,k)
- 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
- qcld(i,k) = qvc(i,k) + zcond(i)
- ENDDO
- !q_sat
- CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:))
- CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:))
-
-
-
! Outputs:
!-------------------------------
Index: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_condensation.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_condensation.F90 (revision 4951)
+++ LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_condensation.F90 (revision 4951)
@@ -0,0 +1,1009 @@
+MODULE lmdz_lscp_condensation
+!----------------------------------------------------------------
+! Module for condensation of clouds routines
+! that are called in LSCP
+
+
+IMPLICIT NONE
+
+CONTAINS
+
+!**********************************************************************************
+SUBROUTINE condensation_lognormal( &
+ klon, temp, qtot, qsat, gamma_cond, ratqs, &
+ keepgoing, cldfra, qincld, qvc)
+
+!----------------------------------------------------------------------
+! This subroutine calculates the formation of clouds, using a
+! statistical cloud scheme. It uses a generalised lognormal distribution
+! of total water in the gridbox
+! See Bony and Emanuel, 2001
+!----------------------------------------------------------------------
+
+USE lmdz_lscp_ini, ONLY: eps
+
+IMPLICIT NONE
+
+!
+! Input
+!
+INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
+!
+REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K]
+REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg]
+REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg]
+REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-]
+REAL, INTENT(IN) , DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-]
+LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed
+!
+! Output
+! NB. those are in INOUT because of the convergence loop on temperature
+! (in some cases, the values are not re-computed) but the values
+! are never used explicitely
+!
+REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-]
+REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg]
+REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg]
+!
+! Local
+!
+INTEGER :: i
+REAL :: pdf_std, pdf_k, pdf_delta
+REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
+!
+!--Loop on klon
+!
+DO i = 1, klon
+
+ !--If a new calculation of the condensation is needed,
+ !--i.e., temperature has not yet converged (or the cloud is
+ !--formed elsewhere)
+ IF (keepgoing(i)) THEN
+
+ pdf_std = ratqs(i) * qtot(i)
+ pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
+ pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
+ pdf_a = pdf_delta / ( pdf_k * SQRT(2.) )
+ pdf_b = pdf_k / (2. * SQRT(2.))
+ pdf_e1 = pdf_a - pdf_b
+ pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
+ pdf_e1 = 1. - ERF(pdf_e1)
+ pdf_e2 = pdf_a + pdf_b
+ pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
+ pdf_e2 = 1. - ERF(pdf_e2)
+
+
+ IF ( pdf_e1 .LT. eps ) THEN
+ cldfra(i) = 0.
+ qincld(i) = qsat(i)
+ !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+ qvc(i) = 0.
+ ELSE
+ cldfra(i) = 0.5 * pdf_e1
+ qincld(i) = qtot(i) * pdf_e2 / pdf_e1
+ !--AB grid-mean vapor in the cloud - we assume saturation adjustment
+ qvc(i) = qsat(i) * cldfra(i)
+ ENDIF
+
+ ENDIF ! end keepgoing
+
+ENDDO ! end loop on i
+END SUBROUTINE condensation_lognormal
+!**********************************************************************************
+
+!**********************************************************************************
+SUBROUTINE condensation_ice_supersat( &
+ klon, dtime, missing_val, pplay, paprsdn, paprsup, &
+ cf_seri, rvc_seri, ratio_qi_qtot, shear, pbl_eps, cell_area, &
+ temp, qtot, qsat, gamma_cond, ratqs, keepgoing, &
+ cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, &
+ dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
+ Tcontr, qcontr, qcontr2, fcontrN, fcontrP, flight_dist, flight_h2o, &
+ dcf_avi, dqi_avi, dqvc_avi)
+
+!----------------------------------------------------------------------
+! This subroutine calculates the formation, evolution and dissipation
+! of clouds, using a process-oriented treatment of the cloud properties
+! (cloud fraction, vapor in the cloud, condensed water in the cloud).
+! It allows for ice supersaturation in cold regions, in clear sky.
+! If ok_unadjusted_clouds, it also allows for sub- and supersaturated
+! cloud water vapors.
+! It also allows for the formation and evolution of condensation trails
+! (contrails) from aviation.
+! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
+! April 2024
+!----------------------------------------------------------------------
+
+USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
+USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
+USE lmdz_lscp_ini, ONLY: eps, temp_nowater, ok_weibull_warm_clouds
+USE lmdz_lscp_ini, ONLY: ok_unadjusted_clouds, iflag_cloud_sublim_pdf
+USE lmdz_lscp_ini, ONLY: lunout
+
+USE lmdz_lscp_ini, ONLY: mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
+ rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, &
+ cond_thresh_pdf_lscp, coef_mixing_lscp, coef_shear_lscp, &
+ chi_mixing_lscp, rho_ice
+
+IMPLICIT NONE
+
+!
+! Input
+!
+INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
+REAL, INTENT(IN) :: dtime ! time step [s]
+REAL, INTENT(IN) :: missing_val ! missing value for output
+!
+REAL, INTENT(IN) , DIMENSION(klon) :: pplay ! layer pressure [Pa]
+REAL, INTENT(IN) , DIMENSION(klon) :: paprsdn ! pressure at the lower interface [Pa]
+REAL, INTENT(IN) , DIMENSION(klon) :: paprsup ! pressure at the upper interface [Pa]
+REAL, INTENT(IN) , DIMENSION(klon) :: cf_seri ! cloud fraction [-]
+REAL, INTENT(IN) , DIMENSION(klon) :: rvc_seri ! gridbox-mean water vapor in cloud [kg/kg]
+REAL, INTENT(IN) , DIMENSION(klon) :: ratio_qi_qtot ! specific ice water content to total specific water ratio [-]
+REAL, INTENT(IN) , DIMENSION(klon) :: shear ! vertical shear [s-1]
+REAL, INTENT(IN) , DIMENSION(klon) :: pbl_eps !
+REAL, INTENT(IN) , DIMENSION(klon) :: cell_area !
+REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K]
+REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg]
+REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg]
+REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-]
+REAL, INTENT(INOUT), DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-]
+LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed
+!
+! Input for aviation
+!
+REAL, INTENT(IN), DIMENSION(klon) :: flight_dist !
+REAL, INTENT(IN), DIMENSION(klon) :: flight_h2o !
+!
+! Output
+! NB. cldfra and qincld should be outputed as cf_seri and qi_seri,
+! or as tendencies (maybe in the future)
+! NB. those are in INOUT because of the convergence loop on temperature
+! (in some cases, the values are not re-computed) but the values
+! are never used explicitely
+!
+REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-]
+REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg]
+REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg]
+REAL, INTENT(INOUT), DIMENSION(klon) :: issrfra ! ISSR fraction [-]
+REAL, INTENT(INOUT), DIMENSION(klon) :: qissr ! gridbox-mean ISSR specific water [kg/kg]
+!
+! Diagnostics for condensation and ice supersaturation
+! NB. idem for the INOUT
+!
+REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_sub ! cloud fraction tendency because of sublimation [s-1]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_con ! cloud fraction tendency because of condensation [s-1]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_mix ! cloud fraction tendency because of cloud mixing [s-1]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_adj ! specific ice content tendency because of temperature adjustment [kg/kg/s]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_sub ! specific ice content tendency because of sublimation [kg/kg/s]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_con ! specific ice content tendency because of condensation [kg/kg/s]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_mix ! specific ice content tendency because of cloud mixing [kg/kg/s]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
+!
+! Diagnostics for aviation
+! NB. idem for the INOUT
+!
+REAL, INTENT(INOUT), DIMENSION(klon) :: Tcontr ! critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) [K]
+REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr !
+REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr2 !
+REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrN !
+REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrP !
+REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_avi ! cloud fraction tendency because of aviation [s-1]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_avi ! specific ice content tendency because of aviation [kg/kg/s]
+REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_avi ! specific cloud water vapor tendency because of aviation [kg/kg/s]
+!
+! Local
+!
+INTEGER :: i
+LOGICAL :: ok_warm_cloud
+REAL, DIMENSION(klon) :: qcld, qzero
+!
+! for lognormal
+REAL :: pdf_std, pdf_k, pdf_delta
+REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
+!
+! for unadjusted clouds
+REAL :: qvapincld, qvapincld_new
+REAL :: qiceincld, qice_ratio
+REAL :: pres_sat, kappa
+REAL :: air_thermal_conduct, water_vapor_diff
+REAL :: r_ice
+!
+! for sublimation
+REAL :: pdf_alpha
+REAL :: dqt_sub
+!
+! for condensation
+REAL, DIMENSION(klon) :: qsatl, dqsatl
+REAL :: clrfra, qclr, sl_clr, rhl_clr
+REAL :: pdf_ratqs, pdf_skew, pdf_scale, pdf_loc
+REAL :: pdf_x, pdf_y, pdf_T
+REAL :: pdf_e3, pdf_e4
+REAL :: cf_cond, qt_cond, dqt_con
+!
+! for mixing
+REAL, DIMENSION(klon) :: subfra, qsub
+REAL :: dqt_mix_sub, dqt_mix_issr
+REAL :: dcf_mix_sub, dcf_mix_issr
+REAL :: dqvc_mix_sub, dqvc_mix_issr
+REAL :: dqt_mix
+REAL :: a_mix, bovera, N_cld_mix, L_mix
+REAL :: envfra_mix, cldfra_mix
+REAL :: L_shear, shear_fra
+REAL :: sigma_mix, issrfra_mix, subfra_mix, qvapinmix
+!
+! for cell properties
+REAL :: rho, rhodz, dz
+!REAL :: V_cell, M_cell
+!
+! for aviation and cell properties
+!REAL :: dqt_avi
+!REAL :: contrail_fra
+!
+!
+!--more local variables for diagnostics
+!--imported from YOMCST.h
+!--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1)
+!--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air
+!--values from Schumann, Meteorol Zeitschrift, 1996
+!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
+!--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
+!REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1)
+!REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1)
+!REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft
+!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
+!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
+!REAL :: Gcontr
+!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2)
+!--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1)
+!REAL :: qsatliqcontr
+
+
+!-----------------------------------------------
+! Initialisations
+!-----------------------------------------------
+
+! Ajout des émissions de H2O dues à l'aviation
+! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
+! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O )
+! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
+! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
+! flight_h2O is in kg H2O / s / cell
+!
+!IF (ok_plane_h2o) THEN
+! q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) )
+!ENDIF
+
+
+qzero(:) = 0.
+
+!--Calculation of qsat w.r.t. liquid
+CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsatl)
+
+!
+!--Loop on klon
+!
+DO i = 1, klon
+
+ !--If a new calculation of the condensation is needed,
+ !--i.e., temperature has not yet converged (or the cloud is
+ !--formed elsewhere)
+ IF (keepgoing(i)) THEN
+
+ !--If the temperature is higher than the threshold below which
+ !--there is no liquid in the gridbox, we activate the usual scheme
+ !--(generalised lognormal from Bony and Emanuel 2001)
+ !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
+ !--all clouds, and the lognormal scheme is not activated
+ IF ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) THEN
+
+ pdf_std = ratqs(i) * qtot(i)
+ pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
+ pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
+ pdf_a = pdf_delta / ( pdf_k * SQRT(2.) )
+ pdf_b = pdf_k / (2. * SQRT(2.))
+ pdf_e1 = pdf_a - pdf_b
+ pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
+ pdf_e1 = 1. - ERF(pdf_e1)
+ pdf_e2 = pdf_a + pdf_b
+ pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
+ pdf_e2 = 1. - ERF(pdf_e2)
+
+
+ IF ( pdf_e1 .LT. eps ) THEN
+ cldfra(i) = 0.
+ qincld(i) = qsat(i)
+ qvc(i) = 0.
+ ELSE
+ cldfra(i) = 0.5 * pdf_e1
+ qincld(i) = qtot(i) * pdf_e2 / pdf_e1
+ qvc(i) = qsat(i) * cldfra(i)
+ ENDIF
+
+ !--If the temperature is lower than temp_nowater, we use the new
+ !--condensation scheme that allows for ice supersaturation
+ ELSE
+
+ !--Initialisation
+ IF ( temp(i) .GT. temp_nowater ) THEN
+ !--If the air mass is warm (liquid water can exist),
+ !--all the memory is lost and the scheme becomes statistical,
+ !--i.e., the sublimation and mixing processes are deactivated,
+ !--and the condensation process is slightly adapted
+ !--This can happen only if ok_weibull_warm_clouds = .TRUE.
+ cldfra(i) = 0.
+ qvc(i) = 0.
+ qcld(i) = 0.
+ ok_warm_cloud = .TRUE.
+ ELSE
+ !--The following barriers ensure that the traced cloud properties
+ !--are consistent. In some rare cases, i.e. the cloud water vapor
+ !--can be greater than the total water in the gridbox
+ cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
+ qcld(i) = MAX(0., MIN(qtot(i), ( ratio_qi_qtot(i) + rvc_seri(i) ) * qtot(i)))
+ qvc(i) = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
+ ok_warm_cloud = .FALSE.
+ ENDIF
+
+ dcf_sub(i) = 0.
+ dqi_sub(i) = 0.
+ dqvc_sub(i) = 0.
+ dqi_adj(i) = 0.
+ dqvc_adj(i) = 0.
+ dcf_con(i) = 0.
+ dqi_con(i) = 0.
+ dqvc_con(i) = 0.
+ dcf_mix(i) = 0.
+ dqi_mix(i) = 0.
+ dqvc_mix(i) = 0.
+
+ issrfra(i) = 0.
+ qissr(i) = 0.
+ subfra(i) = 0.
+ qsub(i) = 0.
+
+ !--Initialisation of the cell properties
+ !--Dry density [kg/m3]
+ rho = pplay(i) / temp(i) / RD
+ !--Dry air mass [kg/m2]
+ rhodz = ( paprsdn(i) - paprsup(i) ) / RG
+ !--Cell thickness [m]
+ dz = rhodz / rho
+ !--Cell volume [m3]
+ !V_cell = dz * cell_area(i)
+ !--Cell dry air mass [kg]
+ !M_cell = rhodz * cell_area(i)
+
+
+ IF ( ok_unadjusted_clouds ) THEN
+ !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
+ !--hypothesis is lost, and the vapor in the cloud is purely prognostic.
+ !
+ !--The deposition equation is
+ !-- dmi/dt = alpha*4pi*C*Svi / ( R*T/esi/Mw/Dv + Ls/ka/T * (Ls*Mw/R/T - 1) )
+ !--from Pruppacher and Klett (2010), where
+ !--mi is the mass of one ice crystal [kg]
+ !--C is the capacitance of an ice crystal [m]
+ !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
+ !--R is the perfect gas constant [J/mol/K]
+ !--T is the temperature [K]
+ !--esi is the saturation pressure w.r.t. ice [Pa]
+ !--Mw is the molar mass of water [kg/mol]
+ !--Dv is the diffusivity of water vapor [m2/s]
+ !--Ls is the specific latent heat of sublimation [J/kg/K]
+ !--ka is the thermal conductivity of dry air [J/m/s/K]
+ !
+ !--We have qice = Nice * mi, where Nice is the ice crystal
+ !--number concentration per kg of moist air
+ !--HYPOTHESIS 1: the ice crystals are spherical, therefore
+ !-- C = r_ice
+ !-- mi = 4/3 * pi * r_ice**3 * rho_ice
+ !--HYPOTHESIS 2: the ice crystals are monodisperse with the
+ !--initial radius r_ice_0.
+ !--NB. this is notably different than the assumption
+ !--of a distributed qice in the cloud made in the sublimation process;
+ !--should it be consistent?
+ !
+ !--As the deposition process does not create new ice crystals,
+ !--and because we assume a same r_ice value for all crystals
+ !--therefore the sublimation process does not destroy ice crystals
+ !--(or, in a limit case, it destroys all ice crystals), then
+ !--Nice is a constant during the sublimation/deposition process.
+ !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI r_ice_0**3 rho_ice )
+ !
+ !--The deposition equation then reads:
+ !-- dqi/dt = alpha*4pi*r_ice*(qvc-qsat)/qsat / ( R*T/esi/Mw/Dv + Ls/ka/T * (Ls*Mw/R/T - 1) ) * Nice
+ !-- dqi/dt = alpha*4pi* (qi / qi_0)**(1/3) *r_ice_0*(qvc-qsat)/qsat &
+ !-- / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
+ !-- * qi_0 / ( 4/3 RPI r_ice_0**3 rho_ice )
+ !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
+ !-- *alpha/qsat/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 r_ice_0**2 rho_ice )
+ !--and we have
+ !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * qi_0**(2/3) / r_ice_0**2
+ !-- dqi/dt = qi**(1/3) * (qvc - qsat) / kappa * qi_0**(2/3) / r_ice_0**2
+ !--where kappa = (2/3*rho_ice)*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
+ !
+ !--This system of equations can be resolved with an exact
+ !--explicit numerical integration, having one variable resolved
+ !--explicitly, the other exactly. The exactly resolved variable is
+ !--the one decreasing, so it is qvc if the process is
+ !--condensation, qi if it is sublimation.
+ !
+ !--kappa is computed as an initialisation constant, as it depends only
+ !--on temperature and other pre-computed values
+ pres_sat = qsat(i) / ( EPS_W + ( 1. - EPS_W ) * qsat(i) ) * pplay(i)
+ !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
+ air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
+ !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
+ water_vapor_diff = 0.211 * ( temp(i) / RTT )**1.94 * ( 101325. / pplay(i) ) * 1.e-4
+ kappa = 2. / 3. * rho_ice * qsat(i) &
+ * ( RV * temp(i) / water_vapor_diff / pres_sat &
+ + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
+ ENDIF
+
+
+ !-------------------------------------------------------------------
+ !-- SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD --
+ !-------------------------------------------------------------------
+
+ !--If there is a cloud
+ IF ( cldfra(i) .GT. eps ) THEN
+
+ qvapincld = qvc(i) / cldfra(i)
+ !--The vapor in cloud cannot be higher than the
+ !--condensation threshold
+ qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i))
+ qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
+
+ !--If the ice water content is too low, the cloud is purely sublimated
+ IF ( qiceincld .LT. eps ) THEN
+ dcf_sub(i) = - cldfra(i)
+ dqvc_sub(i) = - qvc(i)
+ dqi_sub(i) = - ( qcld(i) - qvc(i) )
+
+ cldfra(i) = 0.
+ qcld(i) = 0.
+ qvc(i) = 0.
+
+ !--Else, the cloud is adjusted and sublimated
+ ELSE
+
+ IF ( ok_unadjusted_clouds ) THEN
+ !--Here, the initial vapor in the cloud is qvapincld, and we compute
+ !--the new vapor qvapincld_new
+
+ !--r_ice formula from Sun and Rikus (1999)
+ r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 &
+ + 0.7957 * ( qiceincld * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
+
+ IF ( qvapincld .GE. qsat(i) ) THEN
+ !--If the cloud is initially supersaturated
+ !--Exact explicit integration (qvc exact, qice explicit)
+ qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) &
+ * EXP( - dtime * qiceincld / kappa / r_ice**2. )
+ ELSE
+ !--If the cloud is initially subsaturated
+ !--Exact explicit integration (qice exact, qvc explicit)
+ !--The barrier is set so that the resulting vapor in cloud
+ !--cannot be greater than qsat
+ !--qice_ratio is the ratio between the new ice content and
+ !--the old one, it is comprised between 0 and 1
+ qice_ratio = ( 1. - 2. / 3. / kappa / r_ice**2. * dtime * ( qsat(i) - qvapincld ) )
+
+ IF ( qice_ratio .LT. 0. ) THEN
+ !--If all the ice has been sublimated, we sublimate
+ !--completely the cloud and do not activate the sublimation
+ !--process
+ !--Tendencies and diagnostics
+ dcf_sub(i) = - cldfra(i)
+ dqvc_sub(i) = - qvc(i)
+ dqi_sub(i) = - ( qcld(i) - qvc(i) )
+
+ cldfra(i) = 0.
+ qcld(i) = 0.
+ qvc(i) = 0.
+
+ !--The new vapor in cloud is set to 0 so that the
+ !--sublimation process does not activate
+ qvapincld_new = 0.
+ ELSE
+ !--Else, the sublimation process is activated with the
+ !--diagnosed new cloud water vapor
+ !--The new vapor in the cloud is increased with the
+ !--sublimated ice
+ qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**(3./2.) )
+ !--The new vapor in the cloud cannot be greater than qsat
+ qvapincld_new = MIN(qvapincld_new, qsat(i))
+ ENDIF ! qice_ratio .LT. 0.
+ ENDIF ! qvapincld .GT. qsat(i)
+ ELSE
+ !--We keep the saturation adjustment hypothesis, and the vapor in the
+ !--cloud is set equal to the saturation vapor
+ qvapincld_new = qsat(i)
+ ENDIF ! ok_unadjusted_clouds
+
+ !--If the vapor in cloud is below vapor needed for the cloud to survive
+ IF ( qvapincld .LT. qvapincld_new ) THEN
+ !--Sublimation of the subsaturated cloud
+ !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
+ !--to use.
+ !--iflag = 1 --> uniform distribution
+ !--iflag = 2 --> exponential distribution
+ !--iflag = 3 --> gamma distribution (Karcher et al 2018)
+
+ IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
+ !--Uniform distribution starting at qvapincld
+ pdf_e1 = 1. / ( 2. * qiceincld )
+
+ dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
+ dqt_sub = - cldfra(i) * ( qvapincld_new**2. - qvapincld**2. ) &
+ * pdf_e1 / 2.
+
+ ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
+ !--Exponential distribution starting at qvapincld
+ pdf_alpha = 1. / qiceincld
+ pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )
+
+ dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
+ dqt_sub = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
+ + qvapincld - qvapincld_new * pdf_e1 )
+
+ ELSEIF ( iflag_cloud_sublim_pdf .GE. 3 ) THEN
+ !--Gamma distribution starting at qvapincld
+ pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
+ pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
+ pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
+ pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
+
+ dcf_sub(i) = - cldfra(i) * pdf_e1
+ dqt_sub = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
+ ENDIF
+
+ !--Tendencies and diagnostics
+ dqvc_sub(i) = dcf_sub(i) * qvapincld
+ dqi_sub(i) = dqt_sub - dqvc_sub(i)
+
+ !--Add tendencies
+ cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
+ qcld(i) = MAX(0., qcld(i) + dqt_sub)
+ qvc(i) = MAX(0., qvc(i) + dqvc_sub(i))
+
+ ENDIF ! qvapincld .LT. qvapincld_new
+
+ !--Adjustment of the IWC of the remaining cloud
+ !--to the new vapor in cloud (this can be either positive or negative)
+ dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
+ dqi_adj(i) = - dqvc_adj(i)
+
+ !--Add tendencies
+ !--The vapor in the cloud is updated, but not qcld as it is constant
+ !--through this process, as well as cldfra which is unmodified
+ qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
+
+ ENDIF ! qiceincld .LT. eps
+ ENDIF ! cldfra(i) .GT. eps
+
+
+ !--------------------------------------------------------------------------
+ !-- CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS --
+ !--------------------------------------------------------------------------
+ !--This section relies on a distribution of water in the clear-sky region of
+ !--the mesh.
+
+ !--If there is a clear-sky region
+ IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
+
+ !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
+ qclr = qtot(i) - qcld(i)
+
+ !--New PDF
+ rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
+ rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
+
+ !--Calculation of the properties of the PDF
+ !--Parameterization from IAGOS observations
+ !--pdf_e1 and pdf_e2 will be reused below
+
+ pdf_e1 = beta_pdf_lscp * cell_area(i)**( 1. / 4. ) &
+ * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
+ IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
+ pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
+ ELSE
+ pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
+ ENDIF
+ pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
+ pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
+
+ pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
+ pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
+ pdf_loc = rhl_clr - pdf_scale * pdf_e2
+
+ !--Diagnostics of ratqs
+ ratqs(i) = pdf_std / ( qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. )
+
+ !--Calculation of the newly condensed water and fraction (pronostic)
+ !--Integration of the clear sky PDF between gamma_cond*qsat and +inf
+ !--NB. the calculated values are clear-sky averaged
+
+ pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
+ pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
+ pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
+ pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
+ cf_cond = EXP( - pdf_y )
+ qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
+
+ IF ( ok_warm_cloud ) THEN
+ !--If the statistical scheme is activated, the first
+ !--calculation of the tendencies is kept for the diagnosis of
+ !--the cloud properties, and the condensation process stops
+ cldfra(i) = cf_cond
+ qcld(i) = qt_cond
+ qvc(i) = cldfra(i) * qsat(i)
+
+ ELSEIF ( cf_cond .GT. cond_thresh_pdf_lscp ) THEN
+ !--We check that there is something to condense, if not the
+ !--condensation process stops
+
+ !--Calculation of the limit qclr in cloud value (stored in pdf_e4)
+
+ pdf_e3 = ( LOG( 1. / cond_thresh_pdf_lscp ) ** ( 1. / pdf_alpha ) - pdf_e2 ) &
+ / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. )
+ pdf_e4 = ( 2. * pdf_e3 * pdf_e1 - pdf_x ) &
+ / ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp - 1. )
+ IF ( ( MIN(pdf_e4, rhl_clr) .LT. rhlmid_pdf_lscp ) .OR. ( pdf_e4 .GT. rhl_clr ) ) THEN
+ pdf_e4 = pdf_x / ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp + 1. )
+ ENDIF
+ pdf_e4 = pdf_e4 * qsatl(i) / 100.
+
+ !--Calculation of the final tendencies
+ dcf_con(i) = ( 1. - cldfra(i) ) * ( pdf_e4 - qclr / ( 1. - cldfra(i) ) ) &
+ / ( pdf_e4 - qt_cond / cf_cond )
+ dqt_con = qclr - pdf_e4 * ( 1. - cldfra(i) - dcf_con(i) )
+
+
+ !--Barriers
+ dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i))
+ dqt_con = MIN(dqt_con, qclr)
+
+
+ IF ( ok_unadjusted_clouds ) THEN
+ !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
+ !--the new vapor qvapincld. The timestep is divided by two because we do not
+ !--know when the condensation occurs
+ qvapincld = gamma_cond(i) * qsat(i)
+ qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
+ !--r_ice formula from Sun and Rikus (1999)
+ r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 &
+ + 0.7957 * ( qiceincld * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
+ !--As qvapincld is necessarily greater than qsat, we only
+ !--use the exact explicit formulation
+ !--Exact explicit version
+ qvapincld = qsat(i) + ( qvapincld - qsat(i) ) &
+ * EXP( - dtime / 2. * qiceincld / kappa / r_ice**2. )
+ ELSE
+ !--We keep the saturation adjustment hypothesis, and the vapor in the
+ !--newly formed cloud is set equal to the saturation vapor.
+ qvapincld = qsat(i)
+ ENDIF
+
+ !--Tendency on cloud vapor and diagnostic
+ dqvc_con(i) = qvapincld * dcf_con(i)
+ dqi_con(i) = dqt_con - dqvc_con(i)
+
+ !--Note that the tendencies are NOT added because they are
+ !--added after the mixing process. In the following, the gridbox fraction is
+ !-- 1. - dcf_con(i), and the total water in the gridbox is
+ !-- qtot(i) - dqi_con(i) - dqvc_con(i)
+
+ ENDIF ! ok_warm_cloud, cf_cond .GT. cond_thresh_pdf_lscp
+ ENDIF ! ( 1. - cldfra(i) ) .GT. eps
+
+ !--If there is still clear sky, we diagnose the ISSR
+ !--We recalculte the PDF properties (after the condensation process)
+ IF ( ( ( 1. - dcf_con(i) - cldfra(i) ) .GT. eps ) .AND. .NOT. ok_warm_cloud ) THEN
+ !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
+ qclr = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i)
+
+ !--New PDF
+ rhl_clr = qclr / ( 1. - dcf_con(i) - cldfra(i) ) / qsatl(i) * 100.
+ rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
+
+ !--Calculation of the properties of the PDF
+ !--Parameterization from IAGOS observations
+ !--pdf_e1 and pdf_e2 will be reused below
+
+ pdf_e1 = beta_pdf_lscp * cell_area(i)**( 1. / 4. ) &
+ * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
+ IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
+ pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
+ ELSE
+ pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
+ ENDIF
+ pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
+ pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
+
+ pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
+ pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
+ pdf_loc = rhl_clr - pdf_scale * pdf_e2
+
+ !--We then calculate the part that is greater than qsat
+ !--and consider it supersaturated
+
+ pdf_x = qsat(i) / qsatl(i) * 100.
+ pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
+ pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
+ pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
+ issrfra(i) = EXP( - pdf_y ) * ( 1. - dcf_con(i) - cldfra(i) )
+ qissr(i) = ( pdf_e3 * ( 1. - dcf_con(i) - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
+ ENDIF
+
+ !--Calculation of the subsaturated clear sky fraction and water
+ subfra(i) = 1. - dcf_con(i) - cldfra(i) - issrfra(i)
+ qsub(i) = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i) - qissr(i)
+
+
+ !--------------------------------------
+ !-- CLOUD MIXING --
+ !--------------------------------------
+ !--This process mixes the cloud with its surroundings: the subsaturated clear sky,
+ !--and the supersaturated clear sky. It is activated if the cloud is big enough,
+ !--but does not cover the entire mesh.
+ !
+ IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
+ .AND. .NOT. ok_warm_cloud ) THEN
+
+ !--Initialisation
+ dcf_mix_sub = 0.
+ dqt_mix_sub = 0.
+ dqvc_mix_sub = 0.
+ dcf_mix_issr = 0.
+ dqt_mix_issr = 0.
+ dqvc_mix_issr = 0.
+
+
+ !-- PART 1 - TURBULENT DIFFUSION
+
+ !--Clouds within the mesh are assumed to be ellipses. The length of the
+ !--semi-major axis is a and the length of the semi-minor axis is b.
+ !--N_cld_mix is the number of clouds within the mesh, and
+ !--clouds_perim is the total perimeter of the clouds within the mesh,
+ !--not considering interfaces with other meshes (only the interfaces with clear
+ !--sky are taken into account).
+ !--
+ !--The area of each cloud is A = a * b * RPI,
+ !--and the perimeter of each cloud is
+ !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
+ !--
+ !--With cell_area the area of the cell, we have:
+ !-- cldfra = A * N_cld_mix / cell_area
+ !-- clouds_perim = P * N_cld_mix
+ !--
+ !--We assume that the ratio between b and a is a function of
+ !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
+ !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
+ !--are spherical.
+ !-- b / a = bovera = MAX(0.1, cldfra)
+ bovera = MAX(0.1, cldfra(i))
+ !--The clouds perimeter is imposed using the formula from Morcrette 2012,
+ !--based on observations.
+ !-- clouds_perim_normalized = alpha * cldfra * ( 1. - cldfra ) = clouds_perim / ( norm_constant * cell_area )
+ !--
+ !--With all this, we have
+ !-- a = SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) )
+ !-- P = RPI * a * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
+ !--and therefore, using the perimeter
+ !-- alpha * cldfra * ( 1. - cldfra ) * norm_constant * cell_area
+ !-- = N_cld_mix * RPI &
+ !-- * SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) ) &
+ !-- * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
+ !--and finally
+ N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - dcf_con(i) - cldfra(i) )**2. &
+ * cell_area(i) * ( 1. - dcf_con(i) ) * bovera / RPI &
+ / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2.
+ !--where coef_mix_lscp = ( alpha * norm_constant )**2.
+ !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer
+ !--In particular, it is 0 if cldfra = 1
+ a_mix = SQRT( cell_area(i) * ( 1. - dcf_con(i) ) * cldfra(i) / bovera / N_cld_mix / RPI )
+
+ !--The time required for turbulent diffusion to homogenize a region of size
+ !--L_mix is defined as (L_mix**2./tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
+ !--We compute L_mix and assume that the cloud is mixed over this length
+ L_mix = SQRT( dtime**3. * pbl_eps(i) )
+ !--The mixing length cannot be greater than the semi-minor axis. In this case,
+ !--the entire cloud is mixed.
+ L_mix = MIN(L_mix, a_mix * bovera)
+
+ !--The fraction of clear sky mixed is
+ !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
+ envfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
+ * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2. )
+ !--The fraction of cloudy sky mixed is
+ !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
+ cldfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
+ * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2. )
+
+
+ !-- PART 2 - SHEARING
+
+ !--The clouds are then sheared. We keep the shape and number
+ !--assumptions from before. The clouds are sheared along their
+ !--semi-major axis (a_mix), on the entire cell heigh dz.
+ !--The increase in size is
+ L_shear = coef_shear_lscp * shear(i) * dz * dtime
+ !--therefore, the total increase in fraction is
+ shear_fra = RPI * L_shear * a_mix * bovera * N_cld_mix &
+ / cell_area(i) / ( 1. - dcf_con(i) )
+ !--and the environment and cloud mixed fractions are the same,
+ !--which we add to the previous calculated mixed fractions.
+ !--We therefore assume that the sheared clouds and the turbulent
+ !--mixed clouds are different.
+ envfra_mix = envfra_mix + shear_fra
+ cldfra_mix = cldfra_mix + shear_fra
+
+
+ !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
+
+ !--The environment fraction is allocated to subsaturated sky or supersaturated sky,
+ !--according to the factor sigma_mix. This is computed as the ratio of the
+ !--subsaturated sky fraction to the environment fraction, corrected by a factor
+ !--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
+ !--supersaturated sky is favoured. Physically, this means that it is more likely
+ !--to have supersaturated sky around the cloud than subsaturated sky.
+ sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) )
+ subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) )
+ issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) )
+ cldfra_mix = MIN( cldfra_mix, cldfra(i) )
+
+ !--First, we mix the subsaturated sky (subfra_mix) and the cloud close
+ !--to this fraction (sigma_mix * cldfra_mix).
+ IF ( subfra(i) .GT. eps ) THEN
+
+ IF ( ok_unadjusted_clouds ) THEN
+ !--The subsaturated air is simply added to the cloud,
+ !--with the corresponding cloud fraction
+ !--If the cloud is too subsaturated, the sublimation process
+ !--activated in the following timestep will reduce the cloud
+ !--fraction
+ dcf_mix_sub = subfra_mix
+ dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
+ dqvc_mix_sub = dqt_mix_sub
+
+ ELSE
+ !--We compute the total humidity in the mixed air, which
+ !--can be either sub- or supersaturated.
+ qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
+ + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
+ / ( subfra_mix + cldfra_mix * sigma_mix )
+
+ IF ( qvapinmix .GT. qsat(i) ) THEN
+ !--If the mixed air is supersaturated, we condense the subsaturated
+ !--region which was mixed.
+ dcf_mix_sub = subfra_mix
+ dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
+ dqvc_mix_sub = dcf_mix_sub * qsat(i)
+ ELSE
+ !--Else, we sublimate the cloud which was mixed.
+ dcf_mix_sub = - sigma_mix * cldfra_mix
+ dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
+ dqvc_mix_sub = dcf_mix_sub * qsat(i)
+ ENDIF
+ ENDIF ! ok_unadjusted_clouds
+ ENDIF ! subfra .GT. eps
+
+ !--We then mix the supersaturated sky (issrfra_mix) and the cloud,
+ !--for which the mixed air is always supersatured, therefore
+ !--the cloud necessarily expands
+ IF ( issrfra(i) .GT. eps ) THEN
+
+ IF ( ok_unadjusted_clouds ) THEN
+ !--The ice supersaturated air is simply added to the
+ !--cloud, and supersaturated vapor will be deposited on the
+ !--cloud ice crystals by the deposition process in the
+ !--following timestep
+ dcf_mix_issr = issrfra_mix
+ dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
+ dqvc_mix_issr = dqt_mix_issr
+ ELSE
+ !--In this case, the additionnal vapor condenses
+ dcf_mix_issr = issrfra_mix
+ dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
+ dqvc_mix_issr = dcf_mix_issr * qsat(i)
+ ENDIF ! ok_unadjusted_clouds
+
+
+ ENDIF ! issrfra .GT. eps
+
+ !--Sum up the tendencies from subsaturated sky and supersaturated sky
+ dcf_mix(i) = dcf_mix_sub + dcf_mix_issr
+ dqt_mix = dqt_mix_sub + dqt_mix_issr
+ dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr
+ dqi_mix(i) = dqt_mix - dqvc_mix(i)
+
+ !--Add tendencies
+ issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
+ qissr(i) = MAX(0., qissr(i) - dqt_mix_issr)
+ cldfra(i) = MAX(0., MIN(1. - dcf_con(i), cldfra(i) + dcf_mix(i)))
+ qcld(i) = MAX(0., MIN(qtot(i) - dqi_con(i) - dqvc_con(i), qcld(i) + dqt_mix))
+ qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
+
+ ENDIF ! ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
+
+ !--Finally, we add the tendencies of condensation
+ cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
+ qcld(i) = MIN(qtot(i), qcld(i) + dqvc_con(i) + dqi_con(i))
+ qvc(i) = MIN(qcld(i), qvc(i) + dqvc_con(i))
+
+
+ !----------------------------------------
+ !-- CONTRAILS AND AVIATION --
+ !----------------------------------------
+
+ !--Add a source of cirrus from aviation contrails
+ !IF ( ok_plane_contrail ) THEN
+ ! dcf_avi(i) = 0.
+ ! dqi_avi(i) = 0.
+ ! dqvc_avi(i) = 0.
+ ! ! TODO implement ok_unadjusted_clouds
+ ! IF ( issrfra(i) .GT. eps ) THEN
+ ! contrail_fra = MIN(1., flight_m(i,k) * dtime * contrail_cross_section / V_cell)
+ ! dcf_avi(i) = issrfra(i) * contrail_fra
+ ! dqt_avi = dcf_avi(i) * qissr(i) / issrfra(i)
+ ! dqvc_avi(i) = qsat(i) * dcf_avi(i)
+ !
+ ! !--Add tendencies
+ ! cldfra(i) = cldfra(i) + dcf_avi(i)
+ ! issrfra(i) = issrfra(i) - dcf_avi(i)
+ ! qcld(i) = qcld(i) + dqt_avi
+ ! qvc(i) = qvc(i) + dqvc_avi(i)
+ ! qissr(i) = qissr(i) - dqt_avi
+
+ ! !--Diagnostics
+ ! dqi_avi(i) = dqt_avi - qsat(i) * dcf_avi(i)
+ ! ENDIF
+ ! dcf_avi(i) = dcf_avi(i) / dtime
+ ! dqi_avi(i) = dqi_avi(i) / dtime
+ ! dqvc_avi(i) = dqvc_avi(i) / dtime
+ !ENDIF
+
+
+
+ !-------------------------------------------
+ !-- FINAL BARRIERS AND OUTPUTS --
+ !-------------------------------------------
+
+ IF ( cldfra(i) .LT. eps ) THEN
+ !--If the cloud is too small, it is sublimated.
+ cldfra(i) = 0.
+ qcld(i) = 0.
+ qvc(i) = 0.
+ qincld(i) = qsat(i)
+ ELSE
+ qincld(i) = qcld(i) / cldfra(i)
+ ENDIF ! cldfra .LT. eps
+
+ !--Diagnostics
+ dcf_sub(i) = dcf_sub(i) / dtime
+ dcf_con(i) = dcf_con(i) / dtime
+ dcf_mix(i) = dcf_mix(i) / dtime
+ dqi_adj(i) = dqi_adj(i) / dtime
+ dqi_sub(i) = dqi_sub(i) / dtime
+ dqi_con(i) = dqi_con(i) / dtime
+ dqi_mix(i) = dqi_mix(i) / dtime
+ dqvc_adj(i) = dqvc_adj(i) / dtime
+ dqvc_sub(i) = dqvc_sub(i) / dtime
+ dqvc_con(i) = dqvc_con(i) / dtime
+ dqvc_mix(i) = dqvc_mix(i) / dtime
+
+ ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
+
+ ENDIF ! end keepgoing
+
+ENDDO ! end loop on i
+
+END SUBROUTINE condensation_ice_supersat
+!**********************************************************************************
+
+END MODULE lmdz_lscp_condensation
Index: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_ini.F90 (revision 4951)
@@ -6,6 +6,6 @@
!--------------------
- REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
- !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI)
+ REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W
+ !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W)
REAL, SAVE, PROTECTED :: seuil_neb=0.001 ! cloud fraction threshold: a cloud can precipitate when exceeded
@@ -136,6 +136,68 @@
!$OMP THREADPRIVATE(dist_liq)
- REAL, SAVE, PROTECTED :: tresh_cl=0.0 ! cloud fraction threshold for cloud top search
+ REAL, SAVE, PROTECTED :: tresh_cl=0.0 ! cloud fraction threshold for cloud top search
!$OMP THREADPRIVATE(tresh_cl)
+
+ !--Parameters for condensation and ice supersaturation
+ LOGICAL, SAVE, PROTECTED :: ok_external_lognormal=.FALSE. ! if True, the lognormal condensation scheme is calculated in the lmdz_lscp_condensation routine
+ !$OMP THREADPRIVATE(ok_external_lognormal)
+
+ LOGICAL, SAVE, PROTECTED :: ok_ice_supersat=.FALSE. ! activates the condensation scheme that allows for ice supersaturation
+ !$OMP THREADPRIVATE(ok_ice_supersat)
+
+ LOGICAL, SAVE, PROTECTED :: ok_unadjusted_clouds=.FALSE. ! if True, relax the saturation adjustment assumption for ice clouds
+ !$OMP THREADPRIVATE(ok_unadjusted_clouds)
+
+ LOGICAL, SAVE, PROTECTED :: ok_weibull_warm_clouds=.FALSE. ! if True, the weibull condensation scheme replaces the lognormal condensation scheme at positive temperatures
+ !$OMP THREADPRIVATE(ok_weibull_warm_clouds)
+
+ INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf=3 ! iflag for the distribution of water inside ice clouds
+ !$OMP THREADPRIVATE(iflag_cloud_sublim_pdf)
+
+ REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp=1./3. ! [-] shape factor of the gamma distribution of water inside ice clouds
+ !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
+
+ REAL, SAVE, PROTECTED :: beta_pdf_lscp=2.8e-3 ! []
+ !$OMP THREADPRIVATE(beta_pdf_lscp)
+
+ REAL, SAVE, PROTECTED :: temp_thresh_pdf_lscp=188. ! [K] below this temperature, water vapor is homoegeneously distributed in the clear sky region
+ !$OMP THREADPRIVATE(temp_thresh_pdf_lscp)
+
+ REAL, SAVE, PROTECTED :: rhlmid_pdf_lscp=57.6 ! [%]
+ !$OMP THREADPRIVATE(rhlmid_pdf_lscp)
+
+ REAL, SAVE, PROTECTED :: k0_pdf_lscp=1.75 ! [-] minimum value of the shape parameter for the weibull law
+ !$OMP THREADPRIVATE(k0_pdf_lscp)
+
+ REAL, SAVE, PROTECTED :: kappa_pdf_lscp=0.0147 ! []
+ !$OMP THREADPRIVATE(kappa_pdf_lscp)
+
+ REAL, SAVE, PROTECTED :: rhl0_pdf_lscp=81.5 ! [%]
+ !$OMP THREADPRIVATE(rhl0_pdf_lscp)
+
+ REAL, SAVE, PROTECTED :: cond_thresh_pdf_lscp=1.E-7 ! [-]
+ !$OMP THREADPRIVATE(cond_thresh_pdf_lscp)
+
+ REAL, SAVE, PROTECTED :: a_homofreez=2.349 ! [-]
+ !$OMP THREADPRIVATE(a_homofreez)
+
+ REAL, SAVE, PROTECTED :: b_homofreez=259. ! [K]
+ !$OMP THREADPRIVATE(b_homofreez)
+
+ REAL, SAVE, PROTECTED :: delta_hetfreez=1. ! [-] value between 0 and 1 to simulate for heterogeneous freezing.
+ !$OMP THREADPRIVATE(delta_hetfreez)
+
+ REAL, SAVE, PROTECTED :: coef_mixing_lscp=1e-7 ! [-]
+ !$OMP THREADPRIVATE(coef_mixing_lscp)
+
+ REAL, SAVE, PROTECTED :: coef_shear_lscp=0.1 ! [-]
+ !$OMP THREADPRIVATE(coef_shear_lscp)
+
+ REAL, SAVE, PROTECTED :: chi_mixing_lscp=1.1 ! [-]
+ !$OMP THREADPRIVATE(chi_mixing_lscp)
+
+! REAL, SAVE, PROTECTED :: contrail_cross_section=200000.
+! !$OMP THREADPRIVATE(contrail_cross_section)
+ !--End of the parameters for condensation and ice supersaturation
!--Parameters for poprecip
@@ -225,18 +287,18 @@
CONTAINS
-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, &
- RVTMP2_in, RTT_in,RD_in,RG_in,RPI_in)
+SUBROUTINE lscp_ini(dtime,lunout_in,prt_level_in,ok_ice_supersat_in, iflag_ratqs, fl_cor_ebil_in, &
+ RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RVTMP2_in, &
+ RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in)
USE ioipsl_getin_p_mod, ONLY : getin_p
- USE ice_sursat_mod, ONLY: ice_sursat_init
USE lmdz_cloudth_ini, ONLY : cloudth_ini
REAL, INTENT(IN) :: dtime
INTEGER, INTENT(IN) :: lunout_in,prt_level_in,iflag_ratqs,fl_cor_ebil_in
- LOGICAL, INTENT(IN) :: ok_ice_sursat
+ LOGICAL, INTENT(IN) :: ok_ice_supersat_in
REAL, INTENT(IN) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in
- REAL, INTENT(IN) :: RVTMP2_in, RTT_in, RD_in, RG_in, RPI_in
+ REAL, INTENT(IN) :: RVTMP2_in, RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in
character (len=20) :: modname='lscp_ini_mod'
character (len=80) :: abort_message
@@ -247,6 +309,9 @@
fl_cor_ebil=fl_cor_ebil_in
+ ok_ice_supersat=ok_ice_supersat_in
+
RG=RG_in
RD=RD_in
+ RV=RV_in
RCPD=RCPD_in
RLVTT=RLVTT_in
@@ -257,4 +322,5 @@
RVTMP2=RVTMP2_in
RPI=RPI_in
+ EPS_W=EPS_W_in
@@ -297,4 +363,5 @@
CALL getin_p('tresh_cl',tresh_cl)
CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
+ ! for poprecip
CALL getin_p('ok_poprecip',ok_poprecip)
CALL getin_p('ok_corr_vap_evasub',ok_corr_vap_evasub)
@@ -316,4 +383,24 @@
CALL getin_p('snow_fallspeed_clr',snow_fallspeed_clr)
CALL getin_p('snow_fallspeed_cld',snow_fallspeed_cld)
+ ! for condensation and ice supersaturation
+ CALL getin_p('ok_external_lognormal',ok_external_lognormal)
+ CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds)
+ CALL getin_p('ok_weibull_warm_clouds',ok_weibull_warm_clouds)
+ CALL getin_p('iflag_cloud_sublim_pdf',iflag_cloud_sublim_pdf)
+ CALL getin_p('mu_subl_pdf_lscp',mu_subl_pdf_lscp)
+ CALL getin_p('beta_pdf_lscp',beta_pdf_lscp)
+ CALL getin_p('temp_thresh_pdf_lscp',temp_thresh_pdf_lscp)
+ CALL getin_p('rhlmid_pdf_lscp',rhlmid_pdf_lscp)
+ CALL getin_p('k0_pdf_lscp',k0_pdf_lscp)
+ CALL getin_p('kappa_pdf_lscp',kappa_pdf_lscp)
+ CALL getin_p('rhl0_pdf_lscp',rhl0_pdf_lscp)
+ CALL getin_p('cond_thresh_pdf_lscp',cond_thresh_pdf_lscp)
+ CALL getin_p('a_homofreez',a_homofreez)
+ CALL getin_p('b_homofreez',b_homofreez)
+ CALL getin_p('delta_hetfreez',delta_hetfreez)
+ CALL getin_p('coef_mixing_lscp',coef_mixing_lscp)
+ CALL getin_p('coef_shear_lscp',coef_shear_lscp)
+ CALL getin_p('chi_mixing_lscp',chi_mixing_lscp)
+ !CALL getin_p('contrail_cross_section',contrail_cross_section)
@@ -354,4 +441,5 @@
WRITE(lunout,*) 'lscp_ini, iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp
WRITE(lunout,*) 'lscp_ini, fl_cor_ebil', fl_cor_ebil
+ ! for poprecip
WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
WRITE(lunout,*) 'lscp_ini, ok_corr_vap_evasub', ok_corr_vap_evasub
@@ -367,4 +455,24 @@
WRITE(lunout,*) 'lscp_ini, snow_fallspeed_clr:', snow_fallspeed_clr
WRITE(lunout,*) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld
+ ! for condensation and ice supersaturation
+ WRITE(lunout,*) 'lscp_ini, ok_external_lognormal:', ok_external_lognormal
+ WRITE(lunout,*) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat
+ WRITE(lunout,*) 'lscp_ini, ok_unadjusted_clouds:', ok_unadjusted_clouds
+ WRITE(lunout,*) 'lscp_ini, ok_weibull_warm_clouds:', ok_weibull_warm_clouds
+ WRITE(lunout,*) 'lscp_ini, iflag_cloud_sublim_pdf:', iflag_cloud_sublim_pdf
+ WRITE(lunout,*) 'lscp_ini, mu_subl_pdf_lscp:', mu_subl_pdf_lscp
+ WRITE(lunout,*) 'lscp_ini, beta_pdf_lscp:', beta_pdf_lscp
+ WRITE(lunout,*) 'lscp_ini, temp_thresh_pdf_lscp:', temp_thresh_pdf_lscp
+ WRITE(lunout,*) 'lscp_ini, rhlmid_pdf_lscp:', rhlmid_pdf_lscp
+ WRITE(lunout,*) 'lscp_ini, k0_pdf_lscp:', k0_pdf_lscp
+ WRITE(lunout,*) 'lscp_ini, kappa_pdf_lscp:', kappa_pdf_lscp
+ WRITE(lunout,*) 'lscp_ini, rhl0_pdf_lscp:', rhl0_pdf_lscp
+ WRITE(lunout,*) 'lscp_ini, a_homofreez:', a_homofreez
+ WRITE(lunout,*) 'lscp_ini, b_homofreez:', b_homofreez
+ WRITE(lunout,*) 'lscp_ini, delta_hetfreez', delta_hetfreez
+ WRITE(lunout,*) 'lscp_ini, coef_mixing_lscp:', coef_mixing_lscp
+ WRITE(lunout,*) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp
+ WRITE(lunout,*) 'lscp_ini, chi_mixing_lscp:', chi_mixing_lscp
+! WRITE(lunout,*) 'lscp_ini, contrail_cross_section:', contrail_cross_section
@@ -389,4 +497,21 @@
+ !--Check flags for condensation and ice supersaturation
+ IF ( ok_external_lognormal .AND. ok_ice_supersat ) THEN
+ abort_message = 'in lscp, ok_external_lognormal=y is incompatible with ok_ice_supersat=y'
+ CALL abort_physic (modname,abort_message,1)
+ ENDIF
+
+ IF ( ok_weibull_warm_clouds .AND. .NOT. ok_ice_supersat ) THEN
+ abort_message = 'in lscp, ok_weibull_warm_clouds=y needs ok_ice_supersat=y'
+ CALL abort_physic (modname,abort_message,1)
+ ENDIF
+
+ IF ( ok_unadjusted_clouds .AND. .NOT. ok_ice_supersat ) THEN
+ abort_message = 'in lscp, ok_unadjusted_clouds=y needs ok_ice_supersat=y'
+ CALL abort_physic (modname,abort_message,1)
+ ENDIF
+
+
!AA Temporary initialisation
a_tr_sca(1) = -0.5
@@ -395,6 +520,4 @@
a_tr_sca(4) = -0.5
- IF (ok_ice_sursat) CALL ice_sursat_init()
-
CALL cloudth_ini(iflag_cloudth_vert,iflag_ratqs)
Index: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp_tools.F90 (revision 4951)
@@ -311,5 +311,6 @@
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- use lmdz_lscp_ini, only: iflag_gammasat, t_glace_min, RTT
+ use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
+ use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez
IMPLICIT NONE
@@ -326,7 +327,5 @@
REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi
- REAL fcirrus, fac
- REAL, PARAMETER :: acirrus=2.349
- REAL, PARAMETER :: bcirrus=259.0
+ REAL f_homofreez, fac
INTEGER i
@@ -342,5 +341,5 @@
dgammasatdt(i)=0.
- ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
+ ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. temp_nowater)) THEN
IF (iflag_gammasat .GE. 2) THEN
@@ -355,14 +354,14 @@
IF (iflag_gammasat .GE.1) THEN
- ! homogeneous freezing of aerosols, according to
- ! Koop, 2000 and Karcher 2008, QJRMS
- ! 'Cirrus regime'
- fcirrus=acirrus-temp(i)/bcirrus
- IF (fcirrus .GT. qsl(i)/qsi(i)) THEN
+ ! homogeneous freezing of aerosols, according to
+ ! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
+ ! 'Cirrus regime'
+ f_homofreez=a_homofreez-temp(i)/b_homofreez
+ IF (f_homofreez .GT. qsl(i)/qsi(i)) THEN
gammasat(i)=qsl(i)/qsi(i)
dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
ELSE
- gammasat(i)=fcirrus
- dgammasatdt(i)=-1.0/bcirrus
+ gammasat(i)= (1.-delta_hetfreez) + delta_hetfreez * f_homofreez
+ dgammasatdt(i)=-delta_hetfreez/b_homofreez
ENDIF
@@ -452,4 +451,180 @@
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+FUNCTION GAMMAINC ( p, x )
+
+!*****************************************************************************80
+!
+!! GAMMAINC computes the regularized lower incomplete Gamma Integral
+!
+! Modified:
+!
+! 20 January 2008
+!
+! Author:
+!
+! Original FORTRAN77 version by B Shea.
+! FORTRAN90 version by John Burkardt.
+!
+! Reference:
+!
+! B Shea,
+! Algorithm AS 239:
+! Chi-squared and Incomplete Gamma Integral,
+! Applied Statistics,
+! Volume 37, Number 3, 1988, pages 466-473.
+!
+! Parameters:
+!
+! Input, real X, P, the parameters of the incomplete
+! gamma ratio. 0 <= X, and 0 < P.
+!
+! Output, real GAMMAINC, the value of the incomplete
+! Gamma integral.
+!
+ IMPLICIT NONE
+
+ REAL A
+ REAL AN
+ REAL ARG
+ REAL B
+ REAL C
+ REAL, PARAMETER :: ELIMIT = - 88.0E+00
+ REAL GAMMAINC
+ REAL, PARAMETER :: OFLO = 1.0E+37
+ REAL P
+ REAL, PARAMETER :: PLIMIT = 1000.0E+00
+ REAL PN1
+ REAL PN2
+ REAL PN3
+ REAL PN4
+ REAL PN5
+ REAL PN6
+ REAL RN
+ REAL, PARAMETER :: TOL = 1.0E-14
+ REAL X
+ REAL, PARAMETER :: XBIG = 1.0E+08
+
+ GAMMAINC = 0.0E+00
+
+ IF ( X == 0.0E+00 ) THEN
+ GAMMAINC = 0.0E+00
+ RETURN
+ END IF
+!
+! IF P IS LARGE, USE A NORMAL APPROXIMATION.
+!
+ IF ( PLIMIT < P ) THEN
+
+ PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
+ + 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )
+
+ GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
+ RETURN
+
+ END IF
+!
+! IF X IS LARGE SET GAMMAD = 1.
+!
+ IF ( XBIG < X ) THEN
+ GAMMAINC = 1.0E+00
+ RETURN
+ END IF
+!
+! USE PEARSON'S SERIES EXPANSION.
+! (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
+!
+ IF ( X <= 1.0E+00 .OR. X < P ) THEN
+
+ ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
+ C = 1.0E+00
+ GAMMAINC = 1.0E+00
+ A = P
+
+ DO
+
+ A = A + 1.0E+00
+ C = C * X / A
+ GAMMAINC = GAMMAINC + C
+
+ IF ( C <= TOL ) THEN
+ EXIT
+ END IF
+
+ END DO
+
+ ARG = ARG + LOG ( GAMMAINC )
+
+ IF ( ELIMIT <= ARG ) THEN
+ GAMMAINC = EXP ( ARG )
+ ELSE
+ GAMMAINC = 0.0E+00
+ END IF
+!
+! USE A CONTINUED FRACTION EXPANSION.
+!
+ ELSE
+
+ ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
+ A = 1.0E+00 - P
+ B = A + X + 1.0E+00
+ C = 0.0E+00
+ PN1 = 1.0E+00
+ PN2 = X
+ PN3 = X + 1.0E+00
+ PN4 = X * B
+ GAMMAINC = PN3 / PN4
+
+ DO
+
+ A = A + 1.0E+00
+ B = B + 2.0E+00
+ C = C + 1.0E+00
+ AN = A * C
+ PN5 = B * PN3 - AN * PN1
+ PN6 = B * PN4 - AN * PN2
+
+ IF ( PN6 /= 0.0E+00 ) THEN
+
+ RN = PN5 / PN6
+
+ IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
+ EXIT
+ END IF
+
+ GAMMAINC = RN
+
+ END IF
+
+ PN1 = PN3
+ PN2 = PN4
+ PN3 = PN5
+ PN4 = PN6
+!
+! RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
+!
+ IF ( OFLO <= ABS ( PN5 ) ) THEN
+ PN1 = PN1 / OFLO
+ PN2 = PN2 / OFLO
+ PN3 = PN3 / OFLO
+ PN4 = PN4 / OFLO
+ END IF
+
+ END DO
+
+ ARG = ARG + LOG ( GAMMAINC )
+
+ IF ( ELIMIT <= ARG ) THEN
+ GAMMAINC = 1.0E+00 - EXP ( ARG )
+ ELSE
+ GAMMAINC = 1.0E+00
+ END IF
+
+ END IF
+
+ RETURN
+END FUNCTION GAMMAINC
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
END MODULE lmdz_lscp_tools
Index: LMDZ6/branches/cirrus/libf/phylmd/phyetat0_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/phyetat0_mod.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/phyetat0_mod.F90 (revision 4951)
@@ -21,5 +21,6 @@
du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
- ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, rneb_ancien, radpas, radsol, rain_fall, ratqs, &
+ ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, &
+ cf_ancien, rvc_ancien, radpas, radsol, rain_fall, ratqs, &
rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
@@ -397,5 +398,4 @@
ancien_ok=ancien_ok.AND.phyetat0_get(ql_ancien,"QLANCIEN","QLANCIEN",0.)
ancien_ok=ancien_ok.AND.phyetat0_get(qs_ancien,"QSANCIEN","QSANCIEN",0.)
- ancien_ok=ancien_ok.AND.phyetat0_get(rneb_ancien,"RNEBANCIEN","RNEBANCIEN",0.)
ancien_ok=ancien_ok.AND.phyetat0_get(u_ancien,"UANCIEN","UANCIEN",0.)
ancien_ok=ancien_ok.AND.phyetat0_get(v_ancien,"VANCIEN","VANCIEN",0.)
@@ -412,4 +412,13 @@
prbsw_ancien(:)=0.
ENDIF
+
+ ! cas specifique des variables de la sursaturation par rapport a la glace
+ IF ( ok_ice_supersat ) THEN
+ ancien_ok=ancien_ok.AND.phyetat0_get(cf_ancien,"CFANCIEN","CFANCIEN",0.)
+ ancien_ok=ancien_ok.AND.phyetat0_get(rvc_ancien,"RVCANCIEN","RVCANCIEN",0.)
+ ELSE
+ cf_ancien(:,:)=0.
+ rvc_ancien(:,:)=0.
+ ENDIF
! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
@@ -419,5 +428,4 @@
(maxval(ql_ancien).EQ.minval(ql_ancien)) .OR. &
(maxval(qs_ancien).EQ.minval(qs_ancien)) .OR. &
- (maxval(rneb_ancien).EQ.minval(rneb_ancien)) .OR. &
(maxval(prw_ancien).EQ.minval(prw_ancien)) .OR. &
(maxval(prlw_ancien).EQ.minval(prlw_ancien)) .OR. &
@@ -432,4 +440,11 @@
ancien_ok=.false.
ENDIF
+ ENDIF
+
+ IF ( ok_ice_supersat ) THEN
+ IF ( (maxval(cf_ancien).EQ.minval(cf_ancien)) .OR. &
+ (maxval(rvc_ancien).EQ.minval(rvc_ancien)) ) THEN
+ ancien_ok=.false.
+ ENDIF
ENDIF
Index: LMDZ6/branches/cirrus/libf/phylmd/phyredem.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/phyredem.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/phyredem.F90 (revision 4951)
@@ -19,6 +19,7 @@
zval, rugoro, t_ancien, q_ancien, &
prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
- ql_ancien, qs_ancien, qbs_ancien, rneb_ancien, u_ancien, &
- v_ancien, clwcon, rnebcon, ratqs, pbl_tke, &
+ ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
+ rvc_ancien, u_ancien, v_ancien, &
+ clwcon, rnebcon, ratqs, pbl_tke, &
wake_delta_pbl_tke, zmax0, f0, sig1, w01, &
wake_deltat, wake_deltaq, wake_s, awake_s, &
@@ -252,5 +253,8 @@
ENDIF
- CALL put_field(pass,"RNEBANCIEN", "RNEBANCIEN", rneb_ancien)
+ IF ( ok_ice_supersat ) THEN
+ CALL put_field(pass,"CFANCIEN", "CFANCIEN", cf_ancien)
+ CALL put_field(pass,"RVCANCIEN", "RVCANCIEN", rvc_ancien)
+ ENDIF
CALL put_field(pass,"PRWANCIEN", "PRWANCIEN", prw_ancien)
Index: LMDZ6/branches/cirrus/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/phys_local_var_mod.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/phys_local_var_mod.F90 (revision 4951)
@@ -18,8 +18,6 @@
REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
!$OMP THREADPRIVATE(u_seri, v_seri)
- REAL, SAVE, ALLOCATABLE :: rneb_seri(:,:)
- !$OMP THREADPRIVATE(rneb_seri)
- REAL, SAVE, ALLOCATABLE :: d_rneb_dyn(:,:)
- !$OMP THREADPRIVATE(d_rneb_dyn)
+ REAL, SAVE, ALLOCATABLE :: cf_seri(:,:), rvc_seri(:,:)
+ !$OMP THREADPRIVATE(cf_seri, rvc_seri)
REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),wprime(:,:,:)
!$OMP THREADPRIVATE(l_mixmin, l_mix, wprime)
@@ -38,4 +36,6 @@
REAL, SAVE, ALLOCATABLE :: d_u_dyn(:,:), d_v_dyn(:,:)
!$OMP THREADPRIVATE(d_u_dyn, d_v_dyn)
+ REAL, SAVE, ALLOCATABLE :: d_cf_dyn(:,:), d_rvc_dyn(:,:)
+ !$OMP THREADPRIVATE(d_cf_dyn, d_rvc_dyn)
REAL, SAVE, ALLOCATABLE :: d_tr_dyn(:,:,:)
!$OMP THREADPRIVATE(d_tr_dyn)
@@ -494,44 +494,33 @@
!$OMP THREADPRIVATE(zn2mout)
- REAL, SAVE, ALLOCATABLE :: qclr(:,:)
- !$OMP THREADPRIVATE(qclr)
- REAL, SAVE, ALLOCATABLE :: qcld(:,:)
- !$OMP THREADPRIVATE(qcld)
- REAL, SAVE, ALLOCATABLE :: qss(:,:)
- !$OMP THREADPRIVATE(qss)
- REAL, SAVE, ALLOCATABLE :: qvc(:,:)
- !$OMP THREADPRIVATE(qvc)
- REAL, SAVE, ALLOCATABLE :: rnebclr(:,:)
- !$OMP THREADPRIVATE(rnebclr)
- REAL, SAVE, ALLOCATABLE :: rnebss(:,:)
- !$OMP THREADPRIVATE(rnebss)
- REAL, SAVE, ALLOCATABLE :: gamma_ss(:,:)
- !$OMP THREADPRIVATE(gamma_ss)
- REAL, SAVE, ALLOCATABLE :: N1_ss(:,:)
- !$OMP THREADPRIVATE(N1_ss)
- REAL, SAVE, ALLOCATABLE :: N2_ss(:,:)
- !$OMP THREADPRIVATE(N2_ss)
- REAL, SAVE, ALLOCATABLE :: drneb_sub(:,:)
- !$OMP THREADPRIVATE(drneb_sub)
- REAL, SAVE, ALLOCATABLE :: drneb_con(:,:)
- !$OMP THREADPRIVATE(drneb_con)
- REAL, SAVE, ALLOCATABLE :: drneb_tur(:,:)
- !$OMP THREADPRIVATE(drneb_tur)
- REAL, SAVE, ALLOCATABLE :: drneb_avi(:,:)
- !$OMP THREADPRIVATE(drneb_avi)
- REAL, SAVE, ALLOCATABLE :: zqsatl(:,:)
- !$OMP THREADPRIVATE(zqsatl)
- REAL, SAVE, ALLOCATABLE :: zqsats(:,:)
- !$OMP THREADPRIVATE(zqsats)
- REAL, SAVE, ALLOCATABLE :: Tcontr(:,:)
- !$OMP THREADPRIVATE(Tcontr)
- REAL, SAVE, ALLOCATABLE :: qcontr(:,:)
- !$OMP THREADPRIVATE(qcontr)
- REAL, SAVE, ALLOCATABLE :: qcontr2(:,:)
- !$OMP THREADPRIVATE(qcontr2)
- REAL, SAVE, ALLOCATABLE :: fcontrN(:,:)
- !$OMP THREADPRIVATE(fcontrN)
- REAL, SAVE, ALLOCATABLE :: fcontrP(:,:)
- !$OMP THREADPRIVATE(fcontrP)
+!-- LSCP - condensation and ice supersaturation variables
+ REAL, SAVE, ALLOCATABLE :: qsub(:,:), qissr(:,:), qcld(:,:)
+ !$OMP THREADPRIVATE(qsub, qissr, qcld)
+ REAL, SAVE, ALLOCATABLE :: subfra(:,:), issrfra(:,:)
+ !$OMP THREADPRIVATE(subfra, issrfra)
+ REAL, SAVE, ALLOCATABLE :: gamma_cond(:,:)
+ !$OMP THREADPRIVATE(gamma_cond)
+ REAL, SAVE, ALLOCATABLE :: ratio_qi_qtot(:,:)
+ !$OMP THREADPRIVATE(ratio_qi_qtot)
+ REAL, SAVE, ALLOCATABLE :: dcf_sub(:,:), dcf_con(:,:), dcf_mix(:,:)
+ !$OMP THREADPRIVATE(dcf_sub, dcf_con, dcf_mix)
+ REAL, SAVE, ALLOCATABLE :: dqi_adj(:,:), dqi_sub(:,:), dqi_con(:,:), dqi_mix(:,:)
+ !$OMP THREADPRIVATE(dqi_adj, dqi_sub, dqi_con, dqi_mix)
+ REAL, SAVE, ALLOCATABLE :: dqvc_adj(:,:), dqvc_sub(:,:), dqvc_con(:,:), dqvc_mix(:,:)
+ !$OMP THREADPRIVATE(dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix)
+ REAL, SAVE, ALLOCATABLE :: qsatliq(:,:), qsatice(:,:)
+ !$OMP THREADPRIVATE(qsatliq, qsatice)
+
+!-- LSCP - aviation and contrails variables
+ REAL, SAVE, ALLOCATABLE :: Tcontr(:,:), qcontr(:,:), qcontr2(:,:)
+ !$OMP THREADPRIVATE(Tcontr, qcontr, qcontr2)
+ REAL, SAVE, ALLOCATABLE :: fcontrN(:,:), fcontrP(:,:)
+ !$OMP THREADPRIVATE(fcontrN, fcontrP)
+ REAL, SAVE, ALLOCATABLE :: dcf_avi(:,:), dqi_avi(:,:), dqvc_avi(:,:)
+ !$OMP THREADPRIVATE(dcf_avi, dqi_avi, dqvc_avi)
+ REAL, SAVE, ALLOCATABLE :: flight_dist(:,:), flight_h2o(:,:)
+ !$OMP THREADPRIVATE(flight_dist, flight_h2o)
+
+!-- LSCP - mixed phase clouds variables
REAL, SAVE, ALLOCATABLE :: distcltop(:,:)
!$OMP THREADPRIVATE(distcltop)
@@ -539,6 +528,5 @@
!$OMP THREADPRIVATE(temp_cltop)
-
-!--POPRECIP variables
+!-- LSCP - POPRECIP variables
REAL, SAVE, ALLOCATABLE :: qraindiag(:,:)
!$OMP THREADPRIVATE(qraindiag)
@@ -668,4 +656,5 @@
ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev), qbs_seri(klon,klev))
ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
+ ALLOCATE(cf_seri(klon,klev),rvc_seri(klon,klev))
ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
@@ -678,4 +667,5 @@
ALLOCATE(d_q_dyn2d(klon),d_ql_dyn2d(klon),d_qs_dyn2d(klon), d_qbs_dyn2d(klon))
ALLOCATE(d_u_dyn(klon,klev),d_v_dyn(klon,klev))
+ ALLOCATE(d_cf_dyn(klon,klev),d_rvc_dyn(klon,klev))
ALLOCATE(d_tr_dyn(klon,klev,nbtr)) !RomP
ALLOCATE(d_t_con(klon,klev),d_q_con(klon,klev),d_q_con_zmasse(klon,klev))
@@ -954,15 +944,20 @@
ALLOCATE(zn2mout(klon,6))
-! Supersaturation
- ALLOCATE(rneb_seri(klon,klev))
- ALLOCATE(d_rneb_dyn(klon,klev))
- ALLOCATE(qclr(klon,klev), qcld(klon,klev), qss(klon,klev), qvc(klon,klev))
- ALLOCATE(rnebclr(klon,klev), rnebss(klon,klev), gamma_ss(klon,klev))
- ALLOCATE(N1_ss(klon,klev), N2_ss(klon,klev))
- ALLOCATE(drneb_sub(klon,klev), drneb_con(klon,klev), drneb_tur(klon,klev), drneb_avi(klon,klev))
- ALLOCATE(zqsatl(klon,klev), zqsats(klon,klev))
- ALLOCATE(Tcontr(klon,klev), qcontr(klon,klev), qcontr2(klon,klev), fcontrN(klon,klev), fcontrP(klon,klev))
-
-!--POPRECIP variables
+!-- LSCP - condensation and ice supersaturation variables
+ ALLOCATE(qsub(klon,klev), qissr(klon,klev), qcld(klon,klev))
+ ALLOCATE(subfra(klon,klev), issrfra(klon,klev))
+ ALLOCATE(gamma_cond(klon,klev), ratio_qi_qtot(klon,klev))
+ ALLOCATE(dcf_sub(klon,klev), dcf_con(klon,klev), dcf_mix(klon,klev))
+ ALLOCATE(dqi_adj(klon,klev), dqi_sub(klon,klev), dqi_con(klon,klev), dqi_mix(klon,klev))
+ ALLOCATE(dqvc_adj(klon,klev), dqvc_sub(klon,klev), dqvc_con(klon,klev), dqvc_mix(klon,klev))
+ ALLOCATE(qsatliq(klon,klev), qsatice(klon,klev))
+
+!-- LSCP - aviation and contrails variables
+ ALLOCATE(Tcontr(klon,klev), qcontr(klon,klev), qcontr2(klon,klev))
+ ALLOCATE(fcontrN(klon,klev), fcontrP(klon,klev))
+ ALLOCATE(dcf_avi(klon,klev), dqi_avi(klon,klev), dqvc_avi(klon,klev))
+ ALLOCATE(flight_dist(klon,klev), flight_h2o(klon,klev))
+
+!-- LSCP - POPRECIP variables
ALLOCATE(qraindiag(klon,klev), qsnowdiag(klon,klev))
ALLOCATE(dqreva(klon,klev), dqssub(klon,klev))
@@ -1022,4 +1017,5 @@
DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri)
DEALLOCATE(u_seri,v_seri)
+ DEALLOCATE(cf_seri,rvc_seri)
DEALLOCATE(l_mixmin,l_mix,wprime)
DEALLOCATE(pbl_eps)
@@ -1030,4 +1026,5 @@
DEALLOCATE(d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, d_qbs_dyn2d)
DEALLOCATE(d_u_dyn,d_v_dyn)
+ DEALLOCATE(d_cf_dyn,d_rvc_dyn)
DEALLOCATE(d_tr_dyn) !RomP
DEALLOCATE(d_t_con,d_q_con,d_q_con_zmasse)
@@ -1270,15 +1267,20 @@
DEALLOCATE(zn2mout)
-! Supersaturation
- DEALLOCATE(rneb_seri)
- DEALLOCATE(d_rneb_dyn)
- DEALLOCATE(qclr, qcld, qss, qvc)
- DEALLOCATE(rnebclr, rnebss, gamma_ss)
- DEALLOCATE(N1_ss, N2_ss)
- DEALLOCATE(drneb_sub, drneb_con, drneb_tur, drneb_avi)
- DEALLOCATE(zqsatl, zqsats)
- DEALLOCATE(Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
-
-!--POPRECIP variables
+!-- LSCP - condensation and ice supersaturation variables
+ DEALLOCATE(qsub, qissr, qcld)
+ DEALLOCATE(subfra, issrfra)
+ DEALLOCATE(gamma_cond, ratio_qi_qtot)
+ DEALLOCATE(dcf_sub, dcf_con, dcf_mix)
+ DEALLOCATE(dqi_adj, dqi_sub, dqi_con, dqi_mix)
+ DEALLOCATE(dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix)
+ DEALLOCATE(qsatliq, qsatice)
+
+!-- LSCP - aviation and contrails variables
+ DEALLOCATE(Tcontr, qcontr, qcontr2)
+ DEALLOCATE(fcontrN, fcontrP)
+ DEALLOCATE(dcf_avi, dqi_avi, dqvc_avi)
+ DEALLOCATE(flight_dist, flight_h2o)
+
+!-- LSCP - POPRECIP variables
DEALLOCATE(qraindiag, qsnowdiag)
DEALLOCATE(dqreva, dqssub)
Index: LMDZ6/branches/cirrus/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 4951)
@@ -450,10 +450,8 @@
TYPE(ctrl_out), SAVE :: o_SWdn200clr = ctrl_out((/ 10, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
'SWdn200clr', 'SWdn clear sky at 200mb', 'W/m2', (/ ('', i=1, 10) /))
-
- ! arajouter
- ! type(ctrl_out),save :: o_LWupTOA = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'LWupTOA', &
- ! (/ ('', i=1, 10) /))
- ! type(ctrl_out),save :: o_LWupTOAclr = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'LWupTOAclr', &
- ! (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_LWupTOA = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/), &
+ 'LWupTOA', 'LWup at TOA', 'W/m2', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_LWupTOAclr = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/), &
+ 'LWupTOAclr', 'LWup clear sky at TOA', 'W/m2', (/ ('', i=1, 10) /))
! type(ctrl_out),save :: o_LWdnTOA = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'LWdnTOA', &
! (/ ('', i=1, 10) /))
@@ -2069,51 +2067,73 @@
'albslw3', 'Surface albedo LW3', '-', (/ ('', i=1, 10) /))
-!--aviation & supersaturation
- TYPE(ctrl_out), SAVE :: o_oclr = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'oclr', 'Clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_ocld = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'ocld', 'Cloudy sky total water', 'kg/kg', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_oss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'oss', 'ISSR total water', 'kg/kg', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_ovc = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'ovc', 'In-cloup vapor', 'kg/kg', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_rnebclr = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'rnebclr', 'Clear sky fraction', '-', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_rnebss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'rnebss', 'ISSR fraction', '-', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_rnebseri = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'rnebseri', 'Cloud fraction', '-', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_gammass = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'gammass', 'Gamma supersaturation', '', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_N1_ss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'N1ss', 'N1', '', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_N2_ss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'N2ss', 'N2', '', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_drnebsub = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'drnebsub', 'Cloud fraction change because of sublimation', 's-1', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_drnebcon = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'drnebcon', 'Cloud fraction change because of condensation', 's-1', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_drnebtur = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'drnebtur', 'Cloud fraction change because of turbulence', 's-1', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_drnebavi = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'drnebavi', 'Cloud fraction change because of aviation', 's-1', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_qsatl = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'qsatl', 'Saturation with respect to liquid water', '', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_qsats = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'qsats', 'Saturation with respect to solid water', '', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_flight_m = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'flightm', 'Flown meters', 'm/s/mesh', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
- 'flighth2o', 'H2O flight emission', 'kg H2O/s/mesh', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_Tcontr = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
+!-- LSCP - condensation and ice supersaturation variables
+ TYPE(ctrl_out), SAVE :: o_cfseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'cfseri', 'Cloud fraction', '-', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dcfdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dcfdyn', 'Dynamics cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_rvcseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'rvcseri', 'Cloudy water vapor to total water vapor ratio', '-', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_drvcdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'drvcdyn', 'Dynamics cloudy water vapor to total water vapor ratio tendency', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_qsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'qsub', 'Subsaturated clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_qissr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'qissr', 'Supersaturated clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_qcld = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'qcld', 'Cloudy sky total water', 'kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_subfra = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'subfra', 'Subsaturated clear sky fraction', '-', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_issrfra = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'issrfra', 'Supersaturated clear sky fraction', '-', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_gammacond = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'gammacond', 'Condensation threshold w.r.t. saturation', '-', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dcfsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dcfsub', 'Sublimation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dcfcon = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dcfcon', 'Condensation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dcfmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dcfmix', 'Cloud mixing cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqiadj = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqiadj', 'Temperature adjustment ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqisub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqisub', 'Sublimation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqicon = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqicon', 'Condensation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqimix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqimix', 'Cloud mixing ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqvcadj = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqvcadj', 'Temperature adjustment cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqvcsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqvcsub', 'Sublimation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqvccon = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqvccon', 'Condensation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqvcmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqvcmix', 'Cloud mixing cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_qsatl = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'qsatl', 'Saturation with respect to liquid', 'kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_qsati = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'qsati', 'Saturation with respect to ice', 'kg/kg', (/ ('', i=1, 10) /))
+
+!-- LSCP - aviation variables
+ TYPE(ctrl_out), SAVE :: o_Tcontr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
'Tcontr', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))
- TYPE(ctrl_out), SAVE :: o_qcontr = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
- 'qcontr', 'Specific humidity threshold for contrail formation','Pa', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_qcontr2 = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
- 'qcontr2', 'Specific humidity threshold for contrail formation','kg/kg', (/ ('', i=1, 10) /))
- TYPE(ctrl_out), SAVE :: o_fcontrN = ctrl_out((/ 2, 2, 2, 2, 2, 2, 11, 11, 11, 11/),&
+ TYPE(ctrl_out), SAVE :: o_qcontr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+ 'qcontr', 'Specific humidity threshold for contrail formation', 'Pa', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_qcontr2 = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+ 'qcontr2', 'Specific humidity threshold for contrail formation', 'kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_fcontrN = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
'fcontrN', 'Fraction with non-persistent contrail in clear-sky', '-', (/ ('', i=1,10)/))
- TYPE(ctrl_out), SAVE :: o_fcontrP = ctrl_out((/ 2, 2, 2, 2, 2, 2, 11, 11, 11, 11/),&
+ TYPE(ctrl_out), SAVE :: o_fcontrP = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
'fcontrP', 'Fraction with persistent contrail in ISSR', '-', (/ ('', i=1,10)/))
+ TYPE(ctrl_out), SAVE :: o_dcfavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dcfavi', 'Aviation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqiavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqiavi', 'Aviation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_dqvcavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'dqvcavi', 'Aviation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_flight_dist = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'flightdist', 'Aviation flown distance', 'm/s/mesh', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'flighth2o', 'Aviation H2O flight emission', 'kg H2O/s/mesh', (/ ('', i=1, 10) /))
!!!!!!!!!!!!! Sorties niveaux standards de pression NMC
Index: LMDZ6/branches/cirrus/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/phys_output_write_mod.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/phys_output_write_mod.F90 (revision 4951)
@@ -54,4 +54,5 @@
o_SWdnTOAclr, o_nettop, o_SWup200, &
o_SWup200clr, o_SWdn200, o_SWdn200clr, &
+ o_LWupTOA, o_LWupTOAclr, &
o_LWup200, o_LWup200clr, o_LWdn200, &
o_LWdn200clr, o_sols, o_sols0, &
@@ -216,9 +217,12 @@
o_p_tropopause, o_z_tropopause, o_t_tropopause, &
o_col_O3_strato, o_col_O3_tropo, &
-!--aviation & supersaturation
- o_oclr, o_ocld, o_oss, o_ovc, o_rnebss, o_rnebclr, o_rnebseri, o_gammass, &
- o_N1_ss, o_N2_ss, o_qsatl, o_qsats, &
- o_drnebsub, o_drnebcon, o_drnebtur, o_drnebavi, o_flight_m, o_flight_h2o, &
+!-- LSCP - condensation and ice supersaturation variables
+ o_cfseri, o_dcfdyn, o_rvcseri, o_drvcdyn, &
+ o_qsub, o_qissr, o_qcld, o_subfra, o_issrfra, o_gammacond, &
+ o_dcfsub, o_dcfcon, o_dcfmix, o_dqiadj, o_dqisub, o_dqicon, o_dqimix, &
+ o_dqvcadj, o_dqvcsub, o_dqvccon, o_dqvcmix, o_qsatl, o_qsati, &
+!-- LSCP - aviation variables
o_Tcontr, o_qcontr, o_qcontr2, o_fcontrN, o_fcontrP, &
+ o_dcfavi, o_dqiavi, o_dqvcavi, o_flight_dist, o_flight_h2o, &
!--interactive CO2
o_flx_co2_ocean, o_flx_co2_ocean_cor, &
@@ -258,5 +262,4 @@
#endif
- USE ice_sursat_mod, ONLY: flight_m, flight_h2o
USE lmdz_lscp_ini, ONLY: ok_poprecip
@@ -323,5 +326,4 @@
cldq, flwp, fiwp, ue, ve, uq, vq, &
uwat, vwat, &
- rneb_seri, d_rneb_dyn, &
plcl, plfc, wbeff, convoccur, upwd, dnwd, dnwd0, prw, prlw, prsw, prbsw, water_budget, &
s_pblh, s_pblt, s_lcl, s_therm, uwriteSTD, &
@@ -337,8 +339,12 @@
wdtrainA, wdtrainS, wdtrainM, n2, s2, strig, zcong, zlcl_th, proba_notrig, &
random_notrig, &
- qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
- N1_ss, N2_ss, zqsatl, zqsats, &
+ cf_seri, d_cf_dyn, rvc_seri, d_rvc_dyn, &
+ qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
+ dcf_sub, dcf_con, dcf_mix, &
+ dqi_adj, dqi_sub, dqi_con, dqi_mix, &
+ dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
+ qsatliq, qsatice, &
Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
- drneb_sub, drneb_con, drneb_tur, drneb_avi, &
+ dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
alp_bl_det, alp_bl_fluct_m, alp_bl_conv, &
alp_bl_stat, alp_bl_fluct_tke, slab_wfbils, &
@@ -1066,4 +1072,14 @@
CALL histwrite_phy(o_LWupSFCclr, zx_tmp_fi2d)
CALL histwrite_phy(o_LWdnSFCclr, sollwdownclr)
+
+ IF (vars_defined) THEN
+ zx_tmp_fi2d(:) = lwup(:,klevp1)
+ ENDIF
+ CALL histwrite_phy(o_LWupTOA, zx_tmp_fi2d)
+
+ IF (vars_defined) THEN
+ zx_tmp_fi2d(:) = lwup0(:,klevp1)
+ ENDIF
+ CALL histwrite_phy(o_LWupTOAclr, zx_tmp_fi2d)
IF (type_trac/='inca' .OR. config_inca=='aeNP') THEN
@@ -2006,22 +2022,32 @@
ENDIF
-!--aviation & supersaturation
- IF (ok_ice_sursat) THEN
- CALL histwrite_phy(o_oclr, qclr)
- CALL histwrite_phy(o_ocld, qcld)
- CALL histwrite_phy(o_oss, qss)
- CALL histwrite_phy(o_ovc, qvc)
- CALL histwrite_phy(o_rnebclr, rnebclr)
- CALL histwrite_phy(o_rnebss, rnebss)
- CALL histwrite_phy(o_rnebseri, rneb_seri)
- CALL histwrite_phy(o_gammass, gamma_ss)
- CALL histwrite_phy(o_N1_ss, N1_ss)
- CALL histwrite_phy(o_N2_ss, N2_ss)
- CALL histwrite_phy(o_drnebsub, drneb_sub)
- CALL histwrite_phy(o_drnebcon, drneb_con)
- CALL histwrite_phy(o_drnebtur, drneb_tur)
- CALL histwrite_phy(o_drnebavi, drneb_avi)
- CALL histwrite_phy(o_qsatl, zqsatl)
- CALL histwrite_phy(o_qsats, zqsats)
+!-- LSCP - condensation and supersaturation variables
+ IF (ok_ice_supersat) THEN
+ CALL histwrite_phy(o_cfseri, cf_seri)
+ CALL histwrite_phy(o_dcfdyn, d_cf_dyn)
+ CALL histwrite_phy(o_rvcseri, rvc_seri)
+ CALL histwrite_phy(o_drvcdyn, d_rvc_dyn)
+ CALL histwrite_phy(o_qsub, qsub)
+ CALL histwrite_phy(o_qissr, qissr)
+ CALL histwrite_phy(o_qcld, qcld)
+ CALL histwrite_phy(o_subfra, subfra)
+ CALL histwrite_phy(o_issrfra, issrfra)
+ CALL histwrite_phy(o_gammacond, gamma_cond)
+ CALL histwrite_phy(o_dcfsub, dcf_sub)
+ CALL histwrite_phy(o_dcfcon, dcf_con)
+ CALL histwrite_phy(o_dcfmix, dcf_mix)
+ CALL histwrite_phy(o_dqiadj, dqi_adj)
+ CALL histwrite_phy(o_dqisub, dqi_sub)
+ CALL histwrite_phy(o_dqicon, dqi_con)
+ CALL histwrite_phy(o_dqimix, dqi_mix)
+ CALL histwrite_phy(o_dqvcadj, dqvc_adj)
+ CALL histwrite_phy(o_dqvcsub, dqvc_sub)
+ CALL histwrite_phy(o_dqvccon, dqvc_con)
+ CALL histwrite_phy(o_dqvcmix, dqvc_mix)
+ CALL histwrite_phy(o_qsatl, qsatliq)
+ CALL histwrite_phy(o_qsati, qsatice)
+ ENDIF
+!-- LSCP - aviation variables
+ IF (ok_plane_contrail) THEN
CALL histwrite_phy(o_Tcontr, Tcontr)
CALL histwrite_phy(o_qcontr, qcontr)
@@ -2029,7 +2055,8 @@
CALL histwrite_phy(o_fcontrN, fcontrN)
CALL histwrite_phy(o_fcontrP, fcontrP)
- ENDIF
- IF (ok_plane_contrail) THEN
- CALL histwrite_phy(o_flight_m, flight_m)
+ CALL histwrite_phy(o_dcfavi, dcf_avi)
+ CALL histwrite_phy(o_dqiavi, dqi_avi)
+ CALL histwrite_phy(o_dqvcavi, dqvc_avi)
+ CALL histwrite_phy(o_flight_dist, flight_dist)
ENDIF
IF (ok_plane_h2o) THEN
Index: LMDZ6/branches/cirrus/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/phys_state_var_mod.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/phys_state_var_mod.F90 (revision 4951)
@@ -92,4 +92,6 @@
REAL, ALLOCATABLE, SAVE :: u_ancien(:,:), v_ancien(:,:)
!$OMP THREADPRIVATE(u_ancien, v_ancien)
+ REAL, ALLOCATABLE, SAVE :: cf_ancien(:,:), rvc_ancien(:,:)
+!$OMP THREADPRIVATE(cf_ancien, rvc_ancien)
!!! RomP >>>
REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
@@ -100,6 +102,4 @@
REAL, ALLOCATABLE, SAVE :: clwcon(:,:),rnebcon(:,:)
!$OMP THREADPRIVATE(clwcon,rnebcon)
- REAL, ALLOCATABLE, SAVE :: rneb_ancien(:,:)
-!$OMP THREADPRIVATE(rneb_ancien)
REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:),detrain_cv(:,:),fm_cv(:,:)
!$OMP THREADPRIVATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
@@ -586,9 +586,9 @@
ALLOCATE(prw_ancien(klon), prlw_ancien(klon), prsw_ancien(klon), prbsw_ancien(klon))
ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
+ ALLOCATE(cf_ancien(klon,klev), rvc_ancien(klon,klev))
!!! Rom P >>>
ALLOCATE(tr_ancien(klon,klev,nbtr))
!!! Rom P <<<
ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
- ALLOCATE(rneb_ancien(klon,klev))
ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev),detrain_cv(klon,klev),fm_cv(klon,klev))
ALLOCATE(ratqs(klon,klev))
@@ -809,8 +809,9 @@
DEALLOCATE(zthe, zpic, zval)
DEALLOCATE(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
- DEALLOCATE(qs_ancien, ql_ancien, qbs_ancien, rneb_ancien)
+ DEALLOCATE(qs_ancien, ql_ancien, qbs_ancien)
DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien)
DEALLOCATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
DEALLOCATE(u_ancien, v_ancien)
+ DEALLOCATE(cf_ancien, rvc_ancien)
DEALLOCATE(tr_ancien) !RomP
DEALLOCATE(ratqs, pbl_tke,coefh,coefm)
Index: LMDZ6/branches/cirrus/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/physiq_mod.F90 (revision 4950)
+++ LMDZ6/branches/cirrus/libf/phylmd/physiq_mod.F90 (revision 4951)
@@ -73,5 +73,4 @@
USE tracinca_mod, ONLY: config_inca
USE tropopause_m, ONLY: dyn_tropopause
- USE ice_sursat_mod, ONLY: flight_init, airplane
USE vampir
USE write_field_phy
@@ -157,8 +156,10 @@
! [Variables internes non sauvegardees de la physique]
! Variables locales pour effectuer les appels en serie
- t_seri,q_seri,ql_seri,qs_seri,qbs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
+ t_seri,q_seri,ql_seri,qs_seri,qbs_seri, &
+ u_seri,v_seri,cf_seri,rvc_seri,tr_seri, &
rhcl, &
! Dynamic tendencies (diagnostics)
- 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, &
+ d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn, &
+ d_u_dyn,d_v_dyn,d_cf_dyn,d_rvc_dyn,d_tr_dyn, &
d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d,d_qbs_dyn2d, &
! Physic tendencies
@@ -334,7 +335,12 @@
pfraclr,pfracld, &
distcltop,temp_cltop, &
- zqsatl, zqsats, &
- qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
+ !-- LSCP - condensation and ice supersaturation variables
+ qsub, qissr, qcld, subfra, issrfra, gamma_cond, ratio_qi_qtot, &
+ dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
+ dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
+ !-- LSCP - aviation and contrails variables
Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
+ dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
+ !
cldemi, &
cldfra, cldtau, fiwc, &
@@ -511,6 +517,6 @@
!
! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
- INTEGER,SAVE :: ivap, iliq, isol, irneb, ibs
-!$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
+ INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, irvc
+!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, irvc)
!
!
@@ -1362,6 +1368,7 @@
iliq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
- irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
ibs = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
+ icf = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
+ irvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
! CALL init_etat0_limit_unstruct
! IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
@@ -1416,25 +1423,25 @@
ENDIF
- IF (ok_ice_sursat.AND.(iflag_ice_thermo.EQ.0)) THEN
- WRITE (lunout, *) ' ok_ice_sursat=y requires iflag_ice_thermo=1 as well'
+ IF (ok_ice_supersat.AND.(iflag_ice_thermo.EQ.0)) THEN
+ WRITE (lunout, *) ' ok_ice_supersat=y requires iflag_ice_thermo=1 as well'
abort_message='see above'
CALL abort_physic(modname,abort_message,1)
ENDIF
- IF (ok_ice_sursat.AND.(nqo.LT.4)) THEN
- WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
- '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
+ IF (ok_ice_supersat.AND.(nqo.LT.5)) THEN
+ WRITE (lunout, *) ' ok_ice_supersat=y requires 5 H2O tracers ', &
+ '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c) but nqo=', nqo, '. Might as well stop here.'
abort_message='see above'
CALL abort_physic(modname,abort_message,1)
ENDIF
- IF (ok_plane_h2o.AND..NOT.ok_ice_sursat) THEN
- WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_sursat=y '
+ IF (ok_plane_h2o.AND..NOT.ok_ice_supersat) THEN
+ WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_supersat=y '
abort_message='see above'
CALL abort_physic(modname,abort_message,1)
ENDIF
- IF (ok_plane_contrail.AND..NOT.ok_ice_sursat) THEN
- WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_sursat=y '
+ IF (ok_plane_contrail.AND..NOT.ok_ice_supersat) THEN
+ WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_supersat=y '
abort_message='see above'
CALL abort_physic(modname,abort_message,1)
@@ -1442,5 +1449,5 @@
IF (ok_bs) THEN
- IF ((ok_ice_sursat.AND.nqo .LT.5).OR.(.NOT.ok_ice_sursat.AND.nqo.LT.4)) THEN
+ IF ((ok_ice_supersat.AND.nqo .LT.6).OR.(.NOT.ok_ice_supersat.AND.nqo.LT.4)) THEN
WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
'but nqo=', nqo
@@ -1870,5 +1877,6 @@
& RG,RD,RCPD,RKAPPA,RLVTT,RETV)
CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT)
- CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI)
+ CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_supersat,iflag_ratqs,fl_cor_ebil, &
+ RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RV,RG,RPI,EPS_W)
CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
RVTMP2, RTT,RD,RG, RV, RPI)
@@ -2342,4 +2350,6 @@
sollwdown(:))
+ !--Init for LSCP - condensation
+ ratio_qi_qtot(:,:) = 0.
@@ -2439,16 +2449,17 @@
q_seri(i,k) = qx(i,k,ivap)
ql_seri(i,k) = qx(i,k,iliq)
- qbs_seri(i,k) = 0.
+ qbs_seri(i,k)= 0.
+ cf_seri(i,k) = 0.
+ rvc_seri(i,k)= 0.
!CR: ATTENTION, on rajoute la variable glace
IF (nqo.EQ.2) THEN !--vapour and liquid only
qs_seri(i,k) = 0.
- rneb_seri(i,k) = 0.
ELSE IF (nqo.EQ.3) THEN !--vapour, liquid and ice
qs_seri(i,k) = qx(i,k,isol)
- rneb_seri(i,k) = 0.
- ELSE IF (nqo.GE.4) THEN !--vapour, liquid, ice and rneb and blowing snow
+ ELSE IF (nqo.GE.4) THEN !--vapour, liquid, ice, blowing snow, cloud fraction and cloudy water vapor to total water vapor ratio
qs_seri(i,k) = qx(i,k,isol)
- IF (ok_ice_sursat) THEN
- rneb_seri(i,k) = qx(i,k,irneb)
+ IF (ok_ice_supersat) THEN
+ cf_seri(i,k) = qx(i,k,icf)
+ rvc_seri(i,k) = qx(i,k,irvc)
ENDIF
IF (ok_bs) THEN
@@ -2526,5 +2537,7 @@
d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
- d_qbs_dyn(:,:) = (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
+ d_qbs_dyn(:,:)= (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
+ d_cf_dyn(:,:) = (cf_seri(:,:)-cf_ancien(:,:))/phys_tstep
+ d_rvc_dyn(:,:)= (rvc_seri(:,:)-rvc_ancien(:,:))/phys_tstep
CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
@@ -2538,6 +2551,4 @@
IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
! !! RomP <<<
- !!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
- d_rneb_dyn(:,:)=0.0
ELSE
d_u_dyn(:,:) = 0.0
@@ -2547,4 +2558,7 @@
d_ql_dyn(:,:) = 0.0
d_qs_dyn(:,:) = 0.0
+ d_qbs_dyn(:,:)= 0.0
+ d_cf_dyn(:,:) = 0.0
+ d_rvc_dyn(:,:)= 0.0
d_q_dyn2d(:) = 0.0
d_ql_dyn2d(:) = 0.0
@@ -2554,6 +2568,4 @@
IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
! !! RomP <<<
- d_rneb_dyn(:,:)=0.0
- d_qbs_dyn(:,:)=0.0
ancien_ok = .TRUE.
ENDIF
@@ -2664,4 +2676,20 @@
ENDIF
ENDIF
+
+ !-- Needed for LSCP - condensation and ice supersaturation
+ IF (ok_ice_supersat) THEN
+ DO k = 1, klev
+ DO i = 1, klon
+ IF ( ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) ) .GT. 0. ) THEN
+ ratio_qi_qtot(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
+ rvc_seri(i,k) = rvc_seri(i,k) * q_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
+ ELSE
+ ratio_qi_qtot(i,k) = 0.
+ rvc_seri(i,k) = 0.
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDIF
+
!
! Re-evaporer l'eau liquide nuageuse
@@ -3888,11 +3916,11 @@
!--mise à jour de flight_m et flight_h2o dans leur module
- IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
- CALL airplane(debut,pphis,pplay,paprs,t_seri)
- ENDIF
+ !IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
+ ! CALL airplane(debut,pphis,pplay,paprs,t_seri)
+ !ENDIF
CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
t_seri, q_seri,ptconv,ratqs, &
- d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
+ d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
pfraclr,pfracld, &
radocond, picefra, rain_lsc, snow_lsc, &
@@ -3900,7 +3928,11 @@
prfl, psfl, rhcl, &
zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
- iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, temp_cltop, &
- qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
- Tcontr, qcontr, qcontr2, fcontrN, fcontrP , &
+ iflag_ice_thermo, distcltop, temp_cltop, cell_area, &
+ cf_seri, rvc_seri, u_seri, v_seri, pbl_eps(:,:,is_ave), &
+ qsub, qissr, qcld, subfra, issrfra, gamma_cond, ratio_qi_qtot, &
+ dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
+ dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
+ Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
+ dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
@@ -5623,7 +5655,8 @@
d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
ENDIF
- !--ice_sursat: nqo=4, on ajoute rneb
- IF (nqo.ge.4 .and. ok_ice_sursat) THEN
- d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
+ !--ice_supersat: nqo=5, we add cloud fraction and cloudy water vapor to total water vapor ratio
+ IF (nqo.ge.5 .and. ok_ice_supersat) THEN
+ d_qx(i,k,icf) = ( cf_seri(i,k) - qx(i,k,icf) ) / phys_tstep
+ d_qx(i,k,irvc) = ( rvc_seri(i,k) - qx(i,k,irvc) ) / phys_tstep
ENDIF
@@ -5659,6 +5692,7 @@
ql_ancien(:,:) = ql_seri(:,:)
qs_ancien(:,:) = qs_seri(:,:)
- qbs_ancien(:,:) = qbs_seri(:,:)
- rneb_ancien(:,:) = rneb_seri(:,:)
+ qbs_ancien(:,:)= qbs_seri(:,:)
+ cf_ancien(:,:) = cf_seri(:,:)
+ rvc_ancien(:,:)= rvc_seri(:,:)
CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)