From 664010fc912186865f076e45b37d064ded965a96 Mon Sep 17 00:00:00 2001 From: Gaelle DELAUTIER <gaelle.delautier@meteo.fr> Date: Tue, 15 May 2018 14:26:19 +0200 Subject: [PATCH] S.riette 15/5/2018 new ICE3/ICE4 corresponding to LRED=T --- src/MNH/rain_ice_red.f90 | 1695 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 1695 insertions(+) create mode 100644 src/MNH/rain_ice_red.f90 diff --git a/src/MNH/rain_ice_red.f90 b/src/MNH/rain_ice_red.f90 new file mode 100644 index 000000000..b2d30ddf9 --- /dev/null +++ b/src/MNH/rain_ice_red.f90 @@ -0,0 +1,1695 @@ +!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +! ######spl + MODULE MODI_RAIN_ICE_RED +! ######################## +! +INTERFACE + SUBROUTINE RAIN_ICE_RED ( OSEDIC,HSEDIM, HSUBG_AUCV_RC, OWARM, KKA, KKU, KKL, & + PTSTEP, KRR, LDMICRO, PEXN, & + PDZZ, PRHODJ, PRHODREF, PEXNREF, PPABST, PCIT, PCLDFR,& + PTHT, PRVT, PRCT, PRRT, PRIT, PRST, & + PRGT, PTHS, PRVS, PRCS, PRRS, PRIS, PRSS, PRGS, & + PINPRC,PINPRR, PINPRR3D, PEVAP3D, & + PINPRS, PINPRG, PSIGS, PINDEP, PSEA, PTOWN, & + PRHT, PRHS, PINPRH, PFPR ) +! +! +LOGICAL, INTENT(IN) :: OSEDIC ! Switch for droplet sedim. +CHARACTER(LEN=4), INTENT(IN) :: HSEDIM ! Sedimentation scheme +CHARACTER(LEN=4), INTENT(IN) :: HSUBG_AUCV_RC ! Switch for rc->rr Subgrid autoconversion + ! Kind of Subgrid autoconversion method +LOGICAL, INTENT(IN) :: OWARM ! .TRUE. allows raindrops to + ! form by warm processes + ! (Kessler scheme) +! +INTEGER, INTENT(IN) :: KKA !near ground array index +INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index +INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO +REAL, INTENT(IN) :: PTSTEP ! Double Time step + ! (single if cold start) +INTEGER, INTENT(IN) :: KRR ! Number of moist variable +LOGICAL, DIMENSION(:,:,:), INTENT(IN) :: LDMICRO ! mask to limit computation +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXN ! Exner function +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Layer thikness (m) +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! Dry density * Jacobian +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF! Reference density +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Reference Exner function +REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! absolute pressure at t +! +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PCIT ! Pristine ice n.c. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! Cloud fraction +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Theta at time t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT ! Water vapor m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRRT ! Rain water m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRIT ! Pristine ice m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRST ! Snow/aggregate m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRGT ! Graupel/hail m.r. at t +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PSIGS ! Sigma_s at t +! +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS ! Theta source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRVS ! Water vapor m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRCS ! Cloud water m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRRS ! Rain water m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRIS ! Pristine ice m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRSS ! Snow/aggregate m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRGS ! Graupel m.r. source + +! +REAL, DIMENSION(:,:), INTENT(OUT) :: PINPRC! Cloud instant precip +REAL, DIMENSION(:,:), INTENT(INOUT) :: PINDEP ! Cloud instant deposition +REAL, DIMENSION(:,:), INTENT(OUT) :: PINPRR! Rain instant precip +REAL, DIMENSION(:,:,:),INTENT(OUT) :: PINPRR3D! Rain inst precip 3D +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PEVAP3D! Rain evap profile +REAL, DIMENSION(:,:), INTENT(OUT) :: PINPRS! Snow instant precip +REAL, DIMENSION(:,:), INTENT(OUT) :: PINPRG! Graupel instant precip +REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA ! Sea Mask +REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PTOWN! Fraction that is town +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: PRHS ! Hail m.r. source +REAL, DIMENSION(:,:), OPTIONAL, INTENT(OUT) :: PINPRH! Hail instant precip +REAL, DIMENSION(:,:,:,:), OPTIONAL, INTENT(OUT) :: PFPR ! upper-air precipitation fluxes +! +END SUBROUTINE RAIN_ICE_RED +END INTERFACE +END MODULE MODI_RAIN_ICE_RED +! ######spl + SUBROUTINE RAIN_ICE_RED ( OSEDIC,HSEDIM, HSUBG_AUCV_RC, OWARM, KKA, KKU, KKL, & + PTSTEP, KRR, LDMICRO, PEXN, & + PDZZ, PRHODJ, PRHODREF, PEXNREF, PPABST, PCIT, PCLDFR,& + PTHT, PRVT, PRCT, PRRT, PRIT, PRST, & + PRGT, PTHS, PRVS, PRCS, PRRS, PRIS, PRSS, PRGS, & + PINPRC,PINPRR, PINPRR3D, PEVAP3D, & + PINPRS, PINPRG, PSIGS, PINDEP, PSEA, PTOWN, & + PRHT, PRHS, PINPRH, PFPR ) +! ###################################################################### +! +!!**** * - compute the explicit microphysical sources +!! +!! PURPOSE +!! ------- +!! The purpose of this routine is to compute the slow microphysical sources +!! which can be computed explicitly +!! +!! +!!** METHOD +!! ------ +!! The autoconversion computation follows Kessler (1969). +!! The sedimentation rate is computed with a time spliting technique and +!! an upstream scheme, written as a difference of non-advective fluxes. This +!! source term is added to the future instant ( split-implicit process ). +!! The others microphysical processes are evaluated at the central instant +!! (split-explicit process ): autoconversion, accretion and rain evaporation. +!! These last 3 terms are bounded in order not to create negative values +!! for the water species at the future instant. +!! +!! EXTERNAL +!! -------- +!! None +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_PARAMETERS +!! JPHEXT : Horizontal external points number +!! JPVEXT : Vertical external points number +!! Module MODD_CONF : +!! CCONF configuration of the model for the first time step +!! Module MODD_CST +!! XP00 ! Reference pressure +!! XRD,XRV ! Gaz constant for dry air, vapor +!! XMD,XMV ! Molecular weight for dry air, vapor +!! XCPD ! Cpd (dry air) +!! XCL ! Cl (liquid) +!! XCI ! Ci (solid) +!! XTT ! Triple point temperature +!! XLVTT ! Vaporization heat constant +!! XALPW,XBETAW,XGAMW ! Constants for saturation vapor pressure +!! function over liquid water +!! XALPI,XBETAI,XGAMI ! Constants for saturation vapor pressure +!! function over solid ice +!! Module MODD_BUDGET: +!! NBUMOD : model in which budget is calculated +!! CBUTYPE : type of desired budget +!! 'CART' for cartesian box configuration +!! 'MASK' for budget zone defined by a mask +!! 'NONE' ' for no budget +!! NBUPROCCTR : process counter used for each budget variable +!! LBU_RTH : logical for budget of RTH (potential temperature) +!! .TRUE. = budget of RTH +!! .FALSE. = no budget of RTH +!! LBU_RRV : logical for budget of RRV (water vapor) +!! .TRUE. = budget of RRV +!! .FALSE. = no budget of RRV +!! LBU_RRC : logical for budget of RRC (cloud water) +!! .TRUE. = budget of RRC +!! .FALSE. = no budget of RRC +!! LBU_RRI : logical for budget of RRI (cloud ice) +!! .TRUE. = budget of RRI +!! .FALSE. = no budget of RRI +!! LBU_RRR : logical for budget of RRR (rain water) +!! .TRUE. = budget of RRR +!! .FALSE. = no budget of RRR +!! LBU_RRS : logical for budget of RRS (aggregates) +!! .TRUE. = budget of RRS +!! .FALSE. = no budget of RRS +!! LBU_RRG : logical for budget of RRG (graupeln) +!! .TRUE. = budget of RRG +!! .FALSE. = no budget of RRG +!! +!! REFERENCE +!! --------- +!! +!! Book1 and Book2 of documentation ( routine RAIN_ICE ) +!! +!! AUTHOR +!! ------ +!! J.-P. Pinty * Laboratoire d'Aerologie* +!! +!! MODIFICATIONS +!! ------------- +!! Original 02/11/95 +!! (J.Viviand) 04/02/97 debug accumulated prcipitation & convert +!! precipitation rate in m/s +!! (J.-P. Pinty) 17/02/97 add budget calls +!! (J.-P. Pinty) 17/11/97 set ice sedim. for cirrus ice, reset RCHONI +!! and RRHONG, reverse order for DEALLOCATE +!! (J.-P. Pinty) 11/02/98 correction of the air dynamical viscosity and +!! add advance of the budget calls +!! (J.-P. Pinty) 18/05/98 correction of the air density in the RIAUTS +!! process +!! (J.-P. Pinty) 18/11/98 split the main routine +!! (V. Masson) 18/11/98 bug in IVEC1 and IVEC2 upper limits +!! (J. Escobar & J.-P. Pinty) +!! 11/12/98 contains and rewrite count+pack +!! (J. Stein & J.-P. Pinty) +!! 14/10/99 correction for very small RIT +!! (J. Escobar & J.-P. Pinty) +!! 24/07/00 correction for very samll m.r. in +!! the sedimentation subroutine +!! (M. Tomasini) 11/05/01 Autoconversion of rc into rr modification to take +!! into account the subgrid variance +!! (cf Redelsperger & Sommeria JAS 86) +!! (G. Molinie) 21/05/99 bug in RRCFRIG process, RHODREF**(-1) missing +!! in RSRIMCG +!! (G. Molinie & J.-P. Pinty) +!! 21/06/99 bug in RACCS process +!! (P. Jabouille) 27/05/04 safety test for case where esw/i(T)> pabs (~Z>40km) +!! (J-.P. Chaboureau) 12/02/05 temperature depending ice-to-snow autocon- +! version threshold (Chaboureau and Pinty GRL 2006) +!! (J.-P. Pinty) 01/01/O1 add the hail category and correction of the +!! wet growth rate of the graupeln +!! (S.Remy & C.Lac) 06/06 Add the cloud sedimentation +!! (S.Remy & C.Lac) 06/06 Sedimentation becoming the last process +!! to settle the precipitating species created during the current time step +!! (S.Remy & C.Lac) 06/06 Modification of the algorithm of sedimentation +!! to settle n times the precipitating species created during Dt/n instead +!! of Dt +!! (C.Lac) 11/06 Optimization of the sedimentation loop for NEC +!! (J.Escobar) 18/01/2008 Parallel Bug in Budget when IMICRO >= 1 +!! --> Path inhibit this test by IMICRO >= 0 allway true +!! (Y.Seity) 03/2008 Add Statistic sedimentation +!! (Y.Seity) 10/2009 Added condition for the raindrop accretion of the aggregates +!! into graupeln process (5.2.6) to avoid negative graupel mixing ratio +!! (V.Masson, C.Lac) 09/2010 Correction in split sedimentation for +!! reproducibility +!! (S. Riette) Oct 2010 Better vectorisation of RAIN_ICE_SEDIMENTATION_STAT +!! (Y. Seity), 02-2012 add possibility to run with reversed vertical levels +!! (L. Bengtsson), 02-2013 Passing in land/sea mask and town fraction in +!! order to use different cloud droplet number conc. over +!! land, sea and urban areas in the cloud sedimentation. +!! (D. Degrauwe), 2013-11: Export upper-air precipitation fluxes PFPR. +!! (S. Riette) Nov 2013 Protection against null sigma +!! Juan 24/09/2012: for BUG Pgi rewrite PACK function on mode_pack_pgi +!! (C. Lac) FIT temporal scheme : instant M removed +!! (JP Pinty), 01-2014 : ICE4 : partial reconversion of hail to graupel +!! July, 2015 (O.Nuissier/F.Duffourg) Add microphysics diagnostic for +!! aircraft, ballon and profiler +!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 +!! C.Lac : 10/2016 : add droplet deposition +!! C.Lac : 01/2017 : correction on droplet deposition +!! J.Escobar : 10/2017 : for real*4 , limit exp() in RAIN_ICE_SLOW with XMNH_HUGE_12_LOG +!! (C. Abiven, Y. Léauté, V. Seigner, S. Riette) Phasing of Turner rain subgrid param +!! (S. Riette) Source code split into several files +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_PARAMETERS +USE MODD_CST +USE MODD_RAIN_ICE_DESCR +USE MODD_RAIN_ICE_PARAM +USE MODD_PARAM_ICE +USE MODD_BUDGET +USE MODD_LES +USE MODI_BUDGET +USE MODI_ICE4_RAINFR_VERT +USE MODI_ICE4_SEDIMENTATION_STAT +USE MODI_ICE4_SEDIMENTATION_SPLIT +USE MODI_ICE4_NUCLEATION_WRAPPER +USE MODI_ICE4_TENDENCIES +USE MODE_FMWRIT +USE MODE_ll +USE MODE_MSG +! +#ifdef MNH_PGI +USE MODE_PACK_PGI +#endif +! +IMPLICIT NONE +! +!* 0.1 Declarations of dummy arguments : +! +! +! +LOGICAL, INTENT(IN) :: OSEDIC ! Switch for droplet sedim. +CHARACTER(LEN=4), INTENT(IN) :: HSEDIM ! Sedimentation scheme +CHARACTER(LEN=4), INTENT(IN) :: HSUBG_AUCV_RC ! Kind of Subgrid autoconversion method +LOGICAL, INTENT(IN) :: OWARM ! .TRUE. allows raindrops to + ! form by warm processes + ! (Kessler scheme) +INTEGER, INTENT(IN) :: KKA !near ground array index +INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index +INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO +REAL, INTENT(IN) :: PTSTEP ! Double Time step (single if cold start) +INTEGER, INTENT(IN) :: KRR ! Number of moist variable +LOGICAL, DIMENSION(:,:,:), INTENT(IN) :: LDMICRO ! mask to limit computation +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXN ! Exner function +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Layer thikness (m) +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! Dry density * Jacobian +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF! Reference density +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Reference Exner function +REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! absolute pressure at t +! +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PCIT ! Pristine ice n.c. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! Convective Mass Flux Cloud fraction +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Theta at time t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT ! Water vapor m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRRT ! Rain water m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRIT ! Pristine ice m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRST ! Snow/aggregate m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRGT ! Graupel/hail m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PSIGS ! Sigma_s at t +! +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS ! Theta source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRVS ! Water vapor m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRCS ! Cloud water m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRRS ! Rain water m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRIS ! Pristine ice m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRSS ! Snow/aggregate m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRGS ! Graupel m.r. source +! +REAL, DIMENSION(:,:), INTENT(OUT) :: PINPRC! Cloud instant precip +REAL, DIMENSION(:,:), INTENT(INOUT) :: PINDEP ! Cloud instant deposition +REAL, DIMENSION(:,:), INTENT(OUT) :: PINPRR! Rain instant precip +REAL, DIMENSION(:,:,:),INTENT(OUT) :: PINPRR3D! Rain inst precip 3D +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PEVAP3D! Rain evap profile +REAL, DIMENSION(:,:), INTENT(OUT) :: PINPRS! Snow instant precip +REAL, DIMENSION(:,:), INTENT(OUT) :: PINPRG! Graupel instant precip +REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA ! Sea Mask +REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PTOWN! Fraction that is town +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: PRHS ! Hail m.r. source +REAL, DIMENSION(:,:), OPTIONAL, INTENT(OUT) :: PINPRH! Hail instant precip +REAL, DIMENSION(:,:,:,:), OPTIONAL, INTENT(OUT) :: PFPR ! upper-air precipitation fluxes +! +!* 0.2 Declarations of local variables : +! +INTEGER :: IIB ! Define the domain where is +INTEGER :: IIE ! the microphysical sources have to be computed +INTEGER :: IIT ! +INTEGER :: IJB ! +INTEGER :: IJE ! +INTEGER :: IJT ! +INTEGER :: IKB, IKTB, IKT! +INTEGER :: IKE, IKTE ! +! +!For packing +INTEGER :: IMICRO ! Case r_x>0 locations +INTEGER, DIMENSION(COUNT(LDMICRO)) :: I1,I2,I3 ! Used to replace the COUNT +INTEGER :: JL ! and PACK intrinsics +! +!Arrays for nucleation call outisde of LDMICRO points +REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) :: ZW ! work array +REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) :: ZT ! Temperature +REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)) :: & + & ZZ_RVHENI_MR, & ! heterogeneous nucleation mixing ratio change + & ZZ_RVHENI ! heterogeneous nucleation +REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)) :: ZZ_LVFACT, ZZ_LSFACT +! +!Diagnostics +REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)) :: ZRAINFR, & + & ZHLC_HCF3D,& ! HLCLOUDS cloud fraction in high water content part + & ZHLC_LCF3D,& ! HLCLOUDS cloud fraction in low water content part + & ZHLC_HRC3D,& ! HLCLOUDS cloud water content in high water content + & ZHLC_LRC3D ! HLCLOUDS cloud water content in low water content +REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2)) :: ZINPRI ! Pristine ice instant precip +! +!Packed variables +REAL, DIMENSION(COUNT(LDMICRO)) :: ZRVT, & ! Water vapor m.r. at t + & ZRCT, & ! Cloud water m.r. at t + & ZRRT, & ! Rain water m.r. at t + & ZRIT, & ! Pristine ice m.r. at t + & ZRST, & ! Snow/aggregate m.r. at t + & ZRGT, & ! Graupel m.r. at t + & ZRHT, & ! Hail m.r. at t + & ZCIT, & ! Pristine ice conc. at t + & ZTHT, & ! Potential temperature + & ZRHODREF, & ! RHO Dry REFerence + & ZZT, & ! Temperature + & ZPRES, & ! Pressure + & ZEXN, & ! EXNer Pressure + & ZLSFACT, & ! L_s/(Pi*C_ph) + & ZLVFACT, & ! L_v/(Pi*C_ph) + & ZSIGMA_RC,& ! Standard deviation of rc at time t + & ZCF, & ! Cloud fraction + & ZHLC_HCF, & ! HLCLOUDS : fraction of High Cloud Fraction in grid + & ZHLC_LCF, & ! HLCLOUDS : fraction of Low Cloud Fraction in grid + ! note that ZCF = ZHLC_HCF + ZHLC_LCF + & ZHLC_HRC, & ! HLCLOUDS : LWC that is High LWC in grid + & ZHLC_LRC ! HLCLOUDS : LWC that is Low LWC in grid + ! note that ZRC = ZHLC_HRC + ZHLC_LRC +! +!Output packed tendencies (for budgets only) +REAL, DIMENSION(COUNT(LDMICRO)) :: ZRVHENI_MR, & ! heterogeneous nucleation mixing ratio change + & ZRCHONI, & ! Homogeneous nucleation + & ZRRHONG_MR, & ! Spontaneous freezing mixing ratio change + & ZRVDEPS, & ! Deposition on r_s, + & ZRIAGGS, & ! Aggregation on r_s + & ZRIAUTS, & ! Autoconversion of r_i for r_s production + & ZRVDEPG, & ! Deposition on r_g + & ZRCAUTR, & ! Autoconversion of r_c for r_r production + & ZRCACCR, & ! Accretion of r_c for r_r production + & ZRREVAV, & ! Evaporation of r_r + & ZRIMLTC_MR, & ! Cloud ice melting mixing ratio change + & ZRCBERI, & ! Bergeron-Findeisen effect + & ZRHMLTR, & ! Melting of the hailstones + & ZRSMLTG, & ! Conversion-Melting of the aggregates + & ZRCMLTSR, & ! Cloud droplet collection onto aggregates by positive temperature + & ZRRACCSS, ZRRACCSG, ZRSACCRG, & ! Rain accretion onto the aggregates + & ZRCRIMSS, ZRCRIMSG, ZRSRIMCG, ZRSRIMCG_MR, & ! Cloud droplet riming of the aggregates + & ZRICFRRG, ZRRCFRIG, ZRICFRR, & ! Rain contact freezing + & ZRCWETG, ZRIWETG, ZRRWETG, ZRSWETG, & ! Graupel wet growth + & ZRCDRYG, ZRIDRYG, ZRRDRYG, ZRSDRYG, & ! Graupel dry growth + & ZRWETGH, & ! Conversion of graupel into hail + & ZRWETGH_MR, & ! Conversion of graupel into hail, mr change + & ZRGMLTR, & ! Melting of the graupel + & ZRCWETH, ZRIWETH, ZRSWETH, ZRGWETH, ZRRWETH, & ! Dry growth of hailstone + & ZRCDRYH, ZRIDRYH, ZRSDRYH, ZRRDRYH, ZRGDRYH, & ! Wet growth of hailstone + & ZRDRYHG ! Conversion of hailstone into graupel +! +!Output packed total mixing ratio change (for budgets only) +REAL, DIMENSION(COUNT(LDMICRO)) :: ZTOT_RVHENI, & ! heterogeneous nucleation mixing ratio change + & ZTOT_RCHONI, & ! Homogeneous nucleation + & ZTOT_RRHONG, & ! Spontaneous freezing mixing ratio change + & ZTOT_RVDEPS, & ! Deposition on r_s, + & ZTOT_RIAGGS, & ! Aggregation on r_s + & ZTOT_RIAUTS, & ! Autoconversion of r_i for r_s production + & ZTOT_RVDEPG, & ! Deposition on r_g + & ZTOT_RCAUTR, & ! Autoconversion of r_c for r_r production + & ZTOT_RCACCR, & ! Accretion of r_c for r_r production + & ZTOT_RREVAV, & ! Evaporation of r_r + & ZTOT_RCRIMSS, ZTOT_RCRIMSG, ZTOT_RSRIMCG, & ! Cloud droplet riming of the aggregates + & ZTOT_RIMLTC, & ! Cloud ice melting mixing ratio change + & ZTOT_RCBERI, & ! Bergeron-Findeisen effect + & ZTOT_RHMLTR, & ! Melting of the hailstones + & ZTOT_RSMLTG, & ! Conversion-Melting of the aggregates + & ZTOT_RCMLTSR, & ! Cloud droplet collection onto aggregates by positive temperature + & ZTOT_RRACCSS, ZTOT_RRACCSG, ZTOT_RSACCRG, & ! Rain accretion onto the aggregates + & ZTOT_RICFRRG, ZTOT_RRCFRIG, ZTOT_RICFRR, & ! Rain contact freezing + & ZTOT_RCWETG, ZTOT_RIWETG, ZTOT_RRWETG, ZTOT_RSWETG, & ! Graupel wet growth + & ZTOT_RCDRYG, ZTOT_RIDRYG, ZTOT_RRDRYG, ZTOT_RSDRYG, & ! Graupel dry growth + & ZTOT_RWETGH, & ! Conversion of graupel into hail + & ZTOT_RGMLTR, & ! Melting of the graupel + & ZTOT_RCWETH, ZTOT_RIWETH, ZTOT_RSWETH, ZTOT_RGWETH, ZTOT_RRWETH, & ! Dry growth of hailstone + & ZTOT_RCDRYH, ZTOT_RIDRYH, ZTOT_RSDRYH, ZTOT_RRDRYH, ZTOT_RGDRYH, & ! Wet growth of hailstone + & ZTOT_RDRYHG ! Conversion of hailstone into graupel +! +!For time- or mixing-ratio- splitting +REAL, DIMENSION(COUNT(LDMICRO)) :: Z0RVT, & ! Water vapor m.r. at the beginig of the current loop + & Z0RCT, & ! Cloud water m.r. at the beginig of the current loop + & Z0RRT, & ! Rain water m.r. at the beginig of the current loop + & Z0RIT, & ! Pristine ice m.r. at the beginig of the current loop + & Z0RST, & ! Snow/aggregate m.r. at the beginig of the current loop + & Z0RGT, & ! Graupel m.r. at the beginig of the current loop + & Z0RHT, & ! Hail m.r. at the beginig of the current loop + & ZA_TH, ZA_RV, ZA_RC, ZA_RR, ZA_RI, ZA_RS, ZA_RG, ZA_RH, & + & ZB_TH, ZB_RV, ZB_RC, ZB_RR, ZB_RI, ZB_RS, ZB_RG, ZB_RH +! +!To take into acount external tendencies inside the splitting +REAL, DIMENSION(COUNT(LDMICRO)) :: ZEXT_RV, & ! External tendencie for rv + ZEXT_RC, & ! External tendencie for rc + ZEXT_RR, & ! External tendencie for rr + ZEXT_RI, & ! External tendencie for ri + ZEXT_RS, & ! External tendencie for rs + ZEXT_RG, & ! External tendencie for rg + ZEXT_RH, & ! External tendencie for rh + ZEXT_TH, & ! External tendencie for th + ZEXT_WW ! Working array +LOGICAL :: LEXT_TEND +! +INTEGER, DIMENSION(COUNT(LDMICRO)) :: IITER ! Number of iterations done (with real tendencies computation) +INTEGER :: INB_ITER_MAX ! Maximum number of iterations (with real tendencies computation) +REAL, DIMENSION(COUNT(LDMICRO)) :: ZTIME, & ! Current integration time (starts with 0 and ends with PTSTEP) + & ZMAXTIME, & ! Time on which we can apply the current tendencies + & ZTIME_THRESHOLD, & ! Time to reach threshold + & ZTIME_LASTCALL ! Integration time when last tendecies call has been done +LOGICAL, DIMENSION(COUNT(LDMICRO)) :: LLCOMPUTE ! Points where we must compute tendenceis +LOGICAL :: LSOFT ! Must we really compute tendencies or only adjust them to new T variables +LOGICAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2)):: GDEP +REAL :: ZTSTEP ! length of sub-timestep in case of time splitting +REAL :: ZINV_TSTEP ! Inverse ov PTSTEP +REAL, DIMENSION(COUNT(LDMICRO), 6) :: ZRS_TEND +REAL, DIMENSION(COUNT(LDMICRO), 6) :: ZRG_TEND +REAL, DIMENSION(COUNT(LDMICRO), 8) :: ZRH_TEND +! +!For total tendencies computation +REAL, DIMENSION(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)) :: & + &ZW_RVS, ZW_RCS, ZW_RRS, ZW_RIS, ZW_RSS, ZW_RGS, ZW_RHS, ZW_THS +! +!------------------------------------------------------------------------------- +! +!------------------------------------------------------------------------------- +! +!* 1. COMPUTE THE LOOP BOUNDS +! ----------------------- +! +CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) +IIT=SIZE(PDZZ,1) +IJT=SIZE(PDZZ,2) +IKB=KKA+JPVEXT*KKL +IKE=KKU-JPVEXT*KKL +IKT=SIZE(PDZZ,3) +IKTB=1+JPVEXT +IKTE=IKT-JPVEXT +! +ZINV_TSTEP=1./PTSTEP +LEXT_TEND=.TRUE. +! +ZT(:,:,:) = PTHT(:,:,:) * PEXN(:,:,:) +! LSFACT and LVFACT without exner +IF(KRR==7) THEN + ZZ_LSFACT(:,:,:)=(XLSTT+(XCPV-XCI)*(ZT(:,:,:)-XTT)) & + /( XCPD + XCPV*PRVT(:,:,:) + XCL*(PRCT(:,:,:)+PRRT(:,:,:)) & + + XCI*(PRIT(:,:,:)+PRST(:,:,:)+PRGT(:,:,:)+PRHT(:,:,:))) + ZZ_LVFACT(:,:,:)=(XLVTT+(XCPV-XCL)*(ZT(:,:,:)-XTT)) & + /( XCPD + XCPV*PRVT(:,:,:) + XCL*(PRCT(:,:,:)+PRRT(:,:,:)) & + + XCI*(PRIT(:,:,:)+PRST(:,:,:)+PRGT(:,:,:)+PRHT(:,:,:))) +ELSE + ZZ_LSFACT(:,:,:)=(XLSTT+(XCPV-XCI)*(ZT(:,:,:)-XTT)) & + /( XCPD + XCPV*PRVT(:,:,:) + XCL*(PRCT(:,:,:)+PRRT(:,:,:)) & + + XCI*(PRIT(:,:,:)+PRST(:,:,:)+PRGT(:,:,:))) + ZZ_LVFACT(:,:,:)=(XLVTT+(XCPV-XCL)*(ZT(:,:,:)-XTT)) & + /( XCPD + XCPV*PRVT(:,:,:) + XCL*(PRCT(:,:,:)+PRRT(:,:,:)) & + + XCI*(PRIT(:,:,:)+PRST(:,:,:)+PRGT(:,:,:))) +ENDIF +! +!------------------------------------------------------------------------------- +! +!* 2. COMPUTE THE SEDIMENTATION (RS) SOURCE +! ------------------------------------- +! +IF(.NOT. LSEDIM_AFTER) THEN + ! + !* 2.1 sedimentation + ! + IF(HSEDIM=='STAT') THEN + !SR: It *seems* that we must have two separate calls for ifort + IF(KRR==7) THEN + CALL ICE4_SEDIMENTATION_STAT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKTB, IKTE, IKT, KKL, & + &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, & + &PRHODREF, PPABST, PTHT, PRHODJ, & + &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,& + &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,& + &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, & + &PSEA, PTOWN, & + &PINPRH=PINPRH, PRHT=PRHS*PTSTEP, PRHS=PRHS, PFPR=PFPR) + ELSE + CALL ICE4_SEDIMENTATION_STAT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKTB, IKTE, IKT, KKL, & + &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, & + &PRHODREF, PPABST, PTHT, PRHODJ, & + &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,& + &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,& + &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, & + &PSEA, PTOWN, & + &PFPR=PFPR) + ENDIF + PINPRS(:,:) = PINPRS(:,:) + ZINPRI(:,:) + !No negativity correction here as we apply sedimentation on PR.S*PTSTEP variables + ELSEIF(HSEDIM=='SPLI') THEN + !SR: It *seems* that we must have two separate calls for ifort + IF(KRR==7) THEN + CALL ICE4_SEDIMENTATION_SPLIT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKTB, IKTE, IKT, KKL, & + &PTSTEP, KRR, OSEDIC, LDEPOSC,XVDEPOSC, PDZZ, & + &PRHODREF, PPABST, PTHT, PRHODJ, & + &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,& + &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, & + &PSEA, PTOWN, & + &PINPRH=PINPRH, PRHT=PRHT, PRHS=PRHS, PFPR=PFPR) + ELSE + CALL ICE4_SEDIMENTATION_SPLIT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKTB, IKTE, IKT, KKL, & + &PTSTEP, KRR, OSEDIC, LDEPOSC,XVDEPOSC, PDZZ, & + &PRHODREF, PPABST, PTHT, PRHODJ, & + &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,& + &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, & + &PSEA, PTOWN, & + &PFPR=PFPR) + ENDIF + PINPRS(:,:) = PINPRS(:,:) + ZINPRI(:,:) + !We correct negativities with conservation + !SPLI algorith uses a time-splitting. Inside the loop a temporary m.r. is used. + ! It is initialized with the m.r. at T and is modified by two tendencies: + ! sedimentation tendency and an external tendency which represents all other + ! processes (mainly advection and microphysical processes). If both tendencies + ! are negative, sedimentation can remove a specie at a given sub-timestep. From + ! this point sedimentation stops for the remaining sub-timesteps but the other tendency + ! will be still active and will lead to negative values. + ! We could prevent the algorithm to not consume too much a specie, instead we apply + ! a correction here. + CALL CORRECT_NEGATIVITIES(KRR, PRVS, PRCS, PRRS, & + &PRIS, PRSS, PRGS, & + &PTHS, ZZ_LVFACT, ZZ_LSFACT, PRHS) + ELSE + WRITE(*,*) ' STOP' + WRITE(*,*) ' NO SEDIMENTATION SCHEME FOR HSEDIM=', HSEDIM + CALL PRINT_MSG(NVERB_FATAL,'GEN','RAIN_ICE_RED','') + END IF + ! + !* 2.2 budget storage + ! + IF (LBUDGET_RC .AND. OSEDIC) & + CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:), 7 , 'SEDI_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:), 8 , 'SEDI_BU_RRR') + IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:), 9 , 'SEDI_BU_RRI') + IF (LBUDGET_RS) CALL BUDGET (PRSS(:,:,:)*PRHODJ(:,:,:), 10, 'SEDI_BU_RRS') + IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:), 11, 'SEDI_BU_RRG') + IF ( KRR == 7 .AND. LBUDGET_RH) & + CALL BUDGET (PRHS(:,:,:)*PRHODJ(:,:,:), 12, 'SEDI_BU_RRH') + IF ( LBUDGET_RC .AND. LDEPOSC ) & + CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7 ,'DEPO_BU_RRC') +ENDIF +! +!------------------------------------------------------------------------------- +! +!* 3. PACKING +! -------- +! optimization by looking for locations where +! the microphysical fields are larger than a minimal value only !!! +! +IMICRO=0 +IF(COUNT(LDMICRO)/=0) IMICRO=RAIN_ICE_COUNTJV(LDMICRO(:,:,:), IIT, IJT, IKT, SIZE(I1), I1(:), I2(:), I3(:)) +!Packing +IF(IMICRO>0) THEN + DO JL=1, IMICRO + ZRVT(JL) = PRVT(I1(JL),I2(JL),I3(JL)) + ZRCT(JL) = PRCT(I1(JL),I2(JL),I3(JL)) + ZRRT(JL) = PRRT(I1(JL),I2(JL),I3(JL)) + ZRIT(JL) = PRIT(I1(JL),I2(JL),I3(JL)) + ZRST(JL) = PRST(I1(JL),I2(JL),I3(JL)) + ZRGT(JL) = PRGT(I1(JL),I2(JL),I3(JL)) + ZCIT(JL) = PCIT(I1(JL),I2(JL),I3(JL)) + ZCF(JL) = PCLDFR(I1(JL),I2(JL),I3(JL)) + ZRHODREF(JL) = PRHODREF(I1(JL),I2(JL),I3(JL)) + ZTHT(JL) = PTHT(I1(JL),I2(JL),I3(JL)) + ZPRES(JL) = PPABST(I1(JL),I2(JL),I3(JL)) + ZEXN(JL) = PEXN(I1(JL),I2(JL),I3(JL)) + ENDDO + IF(LEXT_TEND) THEN + DO JL=1, IMICRO + ZEXT_RV(JL) = PRVS(I1(JL),I2(JL),I3(JL)) - ZRVT(JL)*ZINV_TSTEP + ZEXT_RC(JL) = PRCS(I1(JL),I2(JL),I3(JL)) - ZRCT(JL)*ZINV_TSTEP + ZEXT_RR(JL) = PRRS(I1(JL),I2(JL),I3(JL)) - ZRRT(JL)*ZINV_TSTEP + ZEXT_RI(JL) = PRIS(I1(JL),I2(JL),I3(JL)) - ZRIT(JL)*ZINV_TSTEP + ZEXT_RS(JL) = PRSS(I1(JL),I2(JL),I3(JL)) - ZRST(JL)*ZINV_TSTEP + ZEXT_RG(JL) = PRGS(I1(JL),I2(JL),I3(JL)) - ZRGT(JL)*ZINV_TSTEP + ZEXT_TH(JL) = PTHS(I1(JL),I2(JL),I3(JL)) - ZTHT(JL)*ZINV_TSTEP + !The th tendency is not related to a mixing ratio change, there is no exn/exnref issue here + ENDDO + ENDIF + IF(HSUBG_AUCV_RC=='PDF ' .AND. CSUBG_PR_PDF=='SIGM') THEN + DO JL=1, IMICRO + ZSIGMA_RC(JL) = PSIGS(I1(JL),I2(JL),I3(JL))*2. +! ZSIGMA_RC(JL) = MAX(PSIGS(I1(JL),I2(JL),I3(JL)) * 2., 1.E-12) + ENDDO + ENDIF + IF(KRR==7) THEN + DO JL=1, IMICRO + ZRHT(JL) = PRHT(I1(JL),I2(JL),I3(JL)) + ENDDO + IF(LEXT_TEND) THEN + DO JL=1, IMICRO + ZEXT_RH(JL) = PRHS(I1(JL),I2(JL),I3(JL)) - ZRHT(JL)*ZINV_TSTEP + ENDDO + ENDIF + ELSE + ZRHT(:)=0. + IF(LEXT_TEND) ZEXT_RH(:)=0. + ENDIF + IF(LBU_ENABLE) THEN + ZTOT_RVHENI(:)=0. + ZTOT_RCHONI(:)=0. + ZTOT_RRHONG(:)=0. + ZTOT_RVDEPS(:)=0. + ZTOT_RIAGGS(:)=0. + ZTOT_RIAUTS(:)=0. + ZTOT_RVDEPG(:)=0. + ZTOT_RCAUTR(:)=0. + ZTOT_RCACCR(:)=0. + ZTOT_RREVAV(:)=0. + ZTOT_RCRIMSS(:)=0. + ZTOT_RCRIMSG(:)=0. + ZTOT_RSRIMCG(:)=0. + ZTOT_RIMLTC(:)=0. + ZTOT_RCBERI(:)=0. + ZTOT_RHMLTR(:)=0. + ZTOT_RSMLTG(:)=0. + ZTOT_RCMLTSR(:)=0. + ZTOT_RRACCSS(:)=0. + ZTOT_RRACCSG(:)=0. + ZTOT_RSACCRG(:)=0. + ZTOT_RICFRRG(:)=0. + ZTOT_RRCFRIG(:)=0. + ZTOT_RICFRR(:)=0. + ZTOT_RCWETG(:)=0. + ZTOT_RIWETG(:)=0. + ZTOT_RRWETG(:)=0. + ZTOT_RSWETG(:)=0. + ZTOT_RCDRYG(:)=0. + ZTOT_RIDRYG(:)=0. + ZTOT_RRDRYG(:)=0. + ZTOT_RSDRYG(:)=0. + ZTOT_RWETGH(:)=0. + ZTOT_RGMLTR(:)=0. + ZTOT_RCWETH(:)=0. + ZTOT_RIWETH(:)=0. + ZTOT_RSWETH(:)=0. + ZTOT_RGWETH(:)=0. + ZTOT_RRWETH(:)=0. + ZTOT_RCDRYH(:)=0. + ZTOT_RIDRYH(:)=0. + ZTOT_RSDRYH(:)=0. + ZTOT_RRDRYH(:)=0. + ZTOT_RGDRYH(:)=0. + ZTOT_RDRYHG(:)=0. + ENDIF +ENDIF +!------------------------------------------------------------------------------- +! +!* 4. LOOP +! ---- +! +!Maximum number of iterations +!We only count real iterations (those for which we *compute* tendencies) +INB_ITER_MAX=NMAXITER +IF(XTSTEP_TS/=0.)THEN + INB_ITER_MAX=MAX(1, INT(PTSTEP/XTSTEP_TS)) !At least the number of iterations needed for the time-splitting + ZTSTEP=PTSTEP/INB_ITER_MAX + INB_ITER_MAX=MAX(NMAXITER, INB_ITER_MAX) !Fot the case XMRSTEP/=0. at the same time +ENDIF +IITER(:)=0 +ZTIME(:)=0. ! Current integration time (all points may have a different integration time) +DO WHILE(ANY(ZTIME(:)<PTSTEP)) ! Loop to *really* compute tendencies + IF(XMRSTEP/=0.) THEN + ! In this case we need to remember the mixing ratios used to compute the tendencies + ! because when mixing ratio has evolved more than a threshold, we must re-compute tendecies + Z0RVT(:)=ZRVT(:) + Z0RCT(:)=ZRCT(:) + Z0RRT(:)=ZRRT(:) + Z0RIT(:)=ZRIT(:) + Z0RST(:)=ZRST(:) + Z0RGT(:)=ZRGT(:) + Z0RHT(:)=ZRHT(:) + ENDIF + IF(XTSTEP_TS/=0.) THEN + ! In this case we need to remember the time when tendencies were computed + ! because when time has evolved more than a limit, we must re-compute tendecies + ZTIME_LASTCALL(:)=ZTIME(:) + ENDIF + LLCOMPUTE(:)=ZTIME(:)<PTSTEP ! Compuation only for points for which integration time has not reached the timestep + LSOFT=.FALSE. ! We *really* compute the tendencies + WHERE(LLCOMPUTE(:)) + IITER(:)=IITER(:)+1 + END WHERE + DO WHILE(ANY(LLCOMPUTE(:))) ! Loop to adjust tendencies when we cross the 0°C or when a specie disappears + ZZT(:) = ZTHT(:) * ZEXN(:) + IF(KRR==7) THEN + ZLSFACT(:)=(XLSTT+(XCPV-XCI)*(ZZT(:)-XTT)) & + /( (XCPD + XCPV*ZRVT(:) + XCL*(ZRCT(:)+ZRRT(:)) & + + XCI*(ZRIT(:)+ZRST(:)+ZRGT(:)+ZRHT(:)))*ZEXN(:) ) + ZLVFACT(:)=(XLVTT+(XCPV-XCL)*(ZZT(:)-XTT)) & + /( (XCPD + XCPV*ZRVT(:) + XCL*(ZRCT(:)+ZRRT(:)) & + + XCI*(ZRIT(:)+ZRST(:)+ZRGT(:)+ZRHT(:)))*ZEXN(:) ) + ELSE + ZLSFACT(:)=(XLSTT+(XCPV-XCI)*(ZZT(:)-XTT)) & + /( (XCPD + XCPV*ZRVT(:) + XCL*(ZRCT(:)+ZRRT(:)) & + + XCI*(ZRIT(:)+ZRST(:)+ZRGT(:)))*ZEXN(:) ) + ZLVFACT(:)=(XLVTT+(XCPV-XCL)*(ZZT(:)-XTT)) & + /( (XCPD + XCPV*ZRVT(:) + XCL*(ZRCT(:)+ZRRT(:)) & + + XCI*(ZRIT(:)+ZRST(:)+ZRGT(:)))*ZEXN(:) ) + ENDIF + ! + !*** 4.1 Tendecies computation + ! + ! Tendencies are *really* computed when LSOFT==.FALSE. and only adjusted otherwise + CALL ICE4_TENDENCIES(IMICRO, IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKT, KKL, & + &KRR, LSOFT, LLCOMPUTE, & + &OWARM, CSUBG_RC_RR_ACCR, CSUBG_RR_EVAP, HSUBG_AUCV_RC, CSUBG_PR_PDF, & + &ZEXN, ZRHODREF, ZLVFACT, ZLSFACT, LDMICRO, I1, I2, I3, & + &ZPRES, ZCF, ZSIGMA_RC, & + &ZCIT, & + &ZZT, ZTHT, & + &ZRVT, ZRCT, ZRRT, ZRIT, ZRST, ZRGT, ZRHT, PRRT, & + &ZRVHENI_MR, ZRRHONG_MR, ZRIMLTC_MR, ZRSRIMCG_MR, & + &ZRCHONI, ZRVDEPS, ZRIAGGS, ZRIAUTS, ZRVDEPG, & + &ZRCAUTR, ZRCACCR, ZRREVAV, & + &ZRCRIMSS, ZRCRIMSG, ZRSRIMCG, ZRRACCSS, ZRRACCSG, ZRSACCRG, ZRSMLTG, ZRCMLTSR, & + &ZRICFRRG, ZRRCFRIG, ZRICFRR, ZRCWETG, ZRIWETG, ZRRWETG, ZRSWETG, & + &ZRCDRYG, ZRIDRYG, ZRRDRYG, ZRSDRYG, ZRWETGH, ZRWETGH_MR, ZRGMLTR, & + &ZRCWETH, ZRIWETH, ZRSWETH, ZRGWETH, ZRRWETH, & + &ZRCDRYH, ZRIDRYH, ZRSDRYH, ZRRDRYH, ZRGDRYH, ZRDRYHG, ZRHMLTR, & + &ZRCBERI, & + &ZRS_TEND, ZRG_TEND, ZRH_TEND, & + &ZA_TH, ZA_RV, ZA_RC, ZA_RR, ZA_RI, ZA_RS, ZA_RG, ZA_RH, & + &ZB_TH, ZB_RV, ZB_RC, ZB_RR, ZB_RI, ZB_RS, ZB_RG, ZB_RH, & + &ZHLC_HCF, ZHLC_LCF, ZHLC_HRC, ZHLC_LRC, ZRAINFR) + ! External tendencies + IF(LEXT_TEND) THEN + ZA_TH(:) = ZA_TH(:) + ZEXT_TH(:) + ZA_RV(:) = ZA_RV(:) + ZEXT_RV(:) + ZA_RC(:) = ZA_RC(:) + ZEXT_RC(:) + ZA_RR(:) = ZA_RR(:) + ZEXT_RR(:) + ZA_RI(:) = ZA_RI(:) + ZEXT_RI(:) + ZA_RS(:) = ZA_RS(:) + ZEXT_RS(:) + ZA_RG(:) = ZA_RG(:) + ZEXT_RG(:) + ZA_RH(:) = ZA_RH(:) + ZEXT_RH(:) + ENDIF + ! + !*** 4.2 Integration time + ! + ! If we can, we will use these tendecies until the end of the timestep + ZMAXTIME(:)=0. + WHERE(LLCOMPUTE(:)) + ZMAXTIME(:)=PTSTEP-ZTIME(:) ! Remaining time until the end of the timestep + ENDWHERE + + !We need to adjust tendencies when temperature reaches 0 + IF(LFEEDBACKT) THEN + !Is ZB_TH enough to change temperature sign? + WHERE( (ZTHT(:) - XTT/ZEXN(:)) * (ZTHT(:) + ZB_TH(:) - XTT/ZEXN(:)) < 0. ) + ZMAXTIME(:)=0. + ENDWHERE + !Can ZA_TH make temperature change of sign? + ZTIME_THRESHOLD(:)=-1. + WHERE(ABS(ZA_TH(:))>1.E-20) + ZTIME_THRESHOLD(:)=(XTT/ZEXN(:) - ZB_TH(:) - ZTHT(:))/ZA_TH(:) + ENDWHERE + WHERE(ZTIME_THRESHOLD(:)>0.) + ZMAXTIME(:)=MIN(ZMAXTIME(:), ZTIME_THRESHOLD(:)) + ENDWHERE + ENDIF + + !We need to adjust tendencies when a specy disappears + !When a specy is missing, only the external tendencies can be negative (and we must keep track of it) + WHERE(ZA_RV(:)<-1.E-20 .AND. ZRVT(:)>XRTMIN(1)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), -(ZB_RV(:)+ZRVT(:))/ZA_RV(:)) + END WHERE + WHERE(ZA_RC(:)<-1.E-20 .AND. ZRCT(:)>XRTMIN(2)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), -(ZB_RC(:)+ZRCT(:))/ZA_RC(:)) + END WHERE + WHERE(ZA_RR(:)<-1.E-20 .AND. ZRRT(:)>XRTMIN(3)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), -(ZB_RR(:)+ZRRT(:))/ZA_RR(:)) + END WHERE + WHERE(ZA_RI(:)<-1.E-20 .AND. ZRIT(:)>XRTMIN(4)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), -(ZB_RI(:)+ZRIT(:))/ZA_RI(:)) + END WHERE + WHERE(ZA_RS(:)<-1.E-20 .AND. ZRST(:)>XRTMIN(5)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), -(ZB_RS(:)+ZRST(:))/ZA_RS(:)) + END WHERE + WHERE(ZA_RG(:)<-1.E-20 .AND. ZRGT(:)>XRTMIN(6)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), -(ZB_RG(:)+ZRGT(:))/ZA_RG(:)) + END WHERE + IF(KRR==7) THEN + WHERE(ZA_RH(:)<-1.E-20 .AND. ZRHT(:)>XRTMIN(7)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), -(ZB_RH(:)+ZRHT(:))/ZA_RH(:)) + END WHERE + ENDIF + + !We stop when the end of the timestep is reached + WHERE(PTSTEP-ZTIME(:)-ZMAXTIME(:)<=0.) + LLCOMPUTE(:)=.FALSE. + ENDWHERE + + !We must recompute tendencies when the end of the sub-timestep is reached + IF(XTSTEP_TS/=0.) THEN + WHERE(IITER(:)<INB_ITER_MAX .AND. ZTIME(:)+ZMAXTIME(:)>ZTIME_LASTCALL(:)+ZTSTEP) + ZMAXTIME(:)=ZTIME_LASTCALL(:)-ZTIME(:)+ZTSTEP + LLCOMPUTE(:)=.FALSE. + ENDWHERE + ENDIF + + !We must recompute tendencies when the maximum allowed change is reached + !When a specy is missing, only the external tendencies can be active and we do not want to recompute + !the microphysical tendencies when external tendencies are negative (results won't change because specy was already missing) + IF(XMRSTEP/=0.) THEN + ZTIME_THRESHOLD(:)=-1. + WHERE(IITER(:)<INB_ITER_MAX .AND. ABS(ZA_RV(:))>1.E-20) + ZTIME_THRESHOLD(:)=(SIGN(1., ZA_RV(:))*XMRSTEP+Z0RVT(:)-ZRVT(:)-ZB_RV(:))/ZA_RV(:) + ENDWHERE + WHERE(ZTIME_THRESHOLD(:)>=0. .AND. ZTIME_THRESHOLD(:)<ZMAXTIME(:) .AND. & + &(ZRVT(:)>XRTMIN(1) .OR. ZA_RV(:)>0.)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), ZTIME_THRESHOLD(:)) + LLCOMPUTE(:)=.FALSE. + ENDWHERE + + ZTIME_THRESHOLD(:)=-1. + WHERE(IITER(:)<INB_ITER_MAX .AND. ABS(ZA_RC(:))>1.E-20) + ZTIME_THRESHOLD(:)=(SIGN(1., ZA_RC(:))*XMRSTEP+Z0RCT(:)-ZRCT(:)-ZB_RC(:))/ZA_RC(:) + ENDWHERE + WHERE(ZTIME_THRESHOLD(:)>=0. .AND. ZTIME_THRESHOLD(:)<ZMAXTIME(:) .AND. & + &(ZRCT(:)>XRTMIN(2) .OR. ZA_RC(:)>0.)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), ZTIME_THRESHOLD(:)) + LLCOMPUTE(:)=.FALSE. + ENDWHERE + + ZTIME_THRESHOLD(:)=-1. + WHERE(IITER(:)<INB_ITER_MAX .AND. ABS(ZA_RR(:))>1.E-20) + ZTIME_THRESHOLD(:)=(SIGN(1., ZA_RR(:))*XMRSTEP+Z0RRT(:)-ZRRT(:)-ZB_RR(:))/ZA_RR(:) + ENDWHERE + WHERE(ZTIME_THRESHOLD(:)>=0. .AND. ZTIME_THRESHOLD(:)<ZMAXTIME(:) .AND. & + &(ZRRT(:)>XRTMIN(3) .OR. ZA_RR(:)>0.)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), ZTIME_THRESHOLD(:)) + LLCOMPUTE(:)=.FALSE. + ENDWHERE + + ZTIME_THRESHOLD(:)=-1. + WHERE(IITER(:)<INB_ITER_MAX .AND. ABS(ZA_RI(:))>1.E-20) + ZTIME_THRESHOLD(:)=(SIGN(1., ZA_RI(:))*XMRSTEP+Z0RIT(:)-ZRIT(:)-ZB_RI(:))/ZA_RI(:) + ENDWHERE + WHERE(ZTIME_THRESHOLD(:)>=0. .AND. ZTIME_THRESHOLD(:)<ZMAXTIME(:) .AND. & + &(ZRIT(:)>XRTMIN(4) .OR. ZA_RI(:)>0.)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), ZTIME_THRESHOLD(:)) + LLCOMPUTE(:)=.FALSE. + ENDWHERE + + ZTIME_THRESHOLD(:)=-1. + WHERE(IITER(:)<INB_ITER_MAX .AND. ABS(ZA_RS(:))>1.E-20) + ZTIME_THRESHOLD(:)=(SIGN(1., ZA_RS(:))*XMRSTEP+Z0RST(:)-ZRST(:)-ZB_RS(:))/ZA_RS(:) + ENDWHERE + WHERE(ZTIME_THRESHOLD(:)>=0. .AND. ZTIME_THRESHOLD(:)<ZMAXTIME(:) .AND. & + &(ZRST(:)>XRTMIN(5) .OR. ZA_RS(:)>0.)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), ZTIME_THRESHOLD(:)) + LLCOMPUTE(:)=.FALSE. + ENDWHERE + + ZTIME_THRESHOLD(:)=-1. + WHERE(IITER(:)<INB_ITER_MAX .AND. ABS(ZA_RG(:))>1.E-20) + ZTIME_THRESHOLD(:)=(SIGN(1., ZA_RG(:))*XMRSTEP+Z0RGT(:)-ZRGT(:)-ZB_RG(:))/ZA_RG(:) + ENDWHERE + WHERE(ZTIME_THRESHOLD(:)>=0. .AND. ZTIME_THRESHOLD(:)<ZMAXTIME(:) .AND. & + &(ZRGT(:)>XRTMIN(6) .OR. ZA_RG(:)>0.)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), ZTIME_THRESHOLD(:)) + LLCOMPUTE(:)=.FALSE. + ENDWHERE + + IF(KRR==7) THEN + ZTIME_THRESHOLD(:)=-1. + WHERE(IITER(:)<INB_ITER_MAX .AND. ABS(ZA_RH(:))>1.E-20) + ZTIME_THRESHOLD(:)=(SIGN(1., ZA_RH(:))*XMRSTEP+Z0RHT(:)-ZRHT(:)-ZB_RH(:))/ZA_RH(:) + ENDWHERE + WHERE(ZTIME_THRESHOLD(:)>=0. .AND. ZTIME_THRESHOLD(:)<ZMAXTIME(:) .AND. & + &(ZRHT(:)>XRTMIN(7) .OR. ZA_RH(:)>0.)) + ZMAXTIME(:)=MIN(ZMAXTIME(:), ZTIME_THRESHOLD(:)) + LLCOMPUTE(:)=.FALSE. + ENDWHERE + ENDIF + + WHERE(IITER(:)<INB_ITER_MAX .AND. MAX(ABS(ZB_RV(:)), ABS(ZB_RC(:)), ABS(ZB_RR(:)), ABS(ZB_RI(:)), & + ABS(ZB_RS(:)), ABS(ZB_RG(:)), ABS(ZB_RH(:)))>XMRSTEP) + ZMAXTIME(:)=0. + LLCOMPUTE(:)=.FALSE. + ENDWHERE + ENDIF + ! + !*** 4.3 New values of variables for next iteration + ! + ZTHT=ZTHT+ZA_TH(:)*ZMAXTIME(:)+ZB_TH(:) + ZRVT=ZRVT+ZA_RV(:)*ZMAXTIME(:)+ZB_RV(:) + ZRCT=ZRCT+ZA_RC(:)*ZMAXTIME(:)+ZB_RC(:) + ZRRT=ZRRT+ZA_RR(:)*ZMAXTIME(:)+ZB_RR(:) + ZRIT=ZRIT+ZA_RI(:)*ZMAXTIME(:)+ZB_RI(:) + ZRST=ZRST+ZA_RS(:)*ZMAXTIME(:)+ZB_RS(:) + ZRGT=ZRGT+ZA_RG(:)*ZMAXTIME(:)+ZB_RG(:) + IF(KRR==7) ZRHT=ZRHT+ZA_RH(:)*ZMAXTIME(:)+ZB_RH(:) + WHERE(ZRIT(:)==0.) + ZCIT(:) = 0. + END WHERE + ! + !*** 4.4 Mixing ratio change due to each process + ! + IF(LBU_ENABLE) THEN + ZTOT_RVHENI(:)= ZTOT_RVHENI(:) +ZRVHENI_MR(:) + ZTOT_RCHONI(:)= ZTOT_RCHONI(:) +ZRCHONI(:) *ZMAXTIME(:) + ZTOT_RRHONG(:)= ZTOT_RRHONG(:) +ZRRHONG_MR(:) + ZTOT_RVDEPS(:)= ZTOT_RVDEPS(:) +ZRVDEPS(:) *ZMAXTIME(:) + ZTOT_RIAGGS(:)= ZTOT_RIAGGS(:) +ZRIAGGS(:) *ZMAXTIME(:) + ZTOT_RIAUTS(:)= ZTOT_RIAUTS(:) +ZRIAUTS(:) *ZMAXTIME(:) + ZTOT_RVDEPG(:)= ZTOT_RVDEPG(:) +ZRVDEPG(:) *ZMAXTIME(:) + ZTOT_RCAUTR(:)= ZTOT_RCAUTR(:) +ZRCAUTR(:) *ZMAXTIME(:) + ZTOT_RCACCR(:)= ZTOT_RCACCR(:) +ZRCACCR(:) *ZMAXTIME(:) + ZTOT_RREVAV(:)= ZTOT_RREVAV(:) +ZRREVAV(:) *ZMAXTIME(:) + ZTOT_RCRIMSS(:)=ZTOT_RCRIMSS(:)+ZRCRIMSS(:)*ZMAXTIME(:) + ZTOT_RCRIMSG(:)=ZTOT_RCRIMSG(:)+ZRCRIMSG(:)*ZMAXTIME(:) + ZTOT_RSRIMCG(:)=ZTOT_RSRIMCG(:)+ZRSRIMCG(:)*ZMAXTIME(:)+ZRSRIMCG_MR(:) + ZTOT_RRACCSS(:)=ZTOT_RRACCSS(:)+ZRRACCSS(:)*ZMAXTIME(:) + ZTOT_RRACCSG(:)=ZTOT_RRACCSG(:)+ZRRACCSG(:)*ZMAXTIME(:) + ZTOT_RSACCRG(:)=ZTOT_RSACCRG(:)+ZRSACCRG(:)*ZMAXTIME(:) + ZTOT_RSMLTG(:)= ZTOT_RSMLTG(:) +ZRSMLTG(:) *ZMAXTIME(:) + ZTOT_RCMLTSR(:)=ZTOT_RCMLTSR(:)+ZRCMLTSR(:) *ZMAXTIME(:) + ZTOT_RICFRRG(:)=ZTOT_RICFRRG(:)+ZRICFRRG(:)*ZMAXTIME(:) + ZTOT_RRCFRIG(:)=ZTOT_RRCFRIG(:)+ZRRCFRIG(:)*ZMAXTIME(:) + ZTOT_RICFRR(:)= ZTOT_RICFRR(:) +ZRICFRR(:) *ZMAXTIME(:) + ZTOT_RCWETG(:)= ZTOT_RCWETG(:) +ZRCWETG(:) *ZMAXTIME(:) + ZTOT_RIWETG(:)= ZTOT_RIWETG(:) +ZRIWETG(:) *ZMAXTIME(:) + ZTOT_RRWETG(:)= ZTOT_RRWETG(:) +ZRRWETG(:) *ZMAXTIME(:) + ZTOT_RSWETG(:)= ZTOT_RSWETG(:) +ZRSWETG(:) *ZMAXTIME(:) + ZTOT_RWETGH(:)= ZTOT_RWETGH(:) +ZRWETGH(:) *ZMAXTIME(:)+ZRWETGH_MR(:) + ZTOT_RCDRYG(:)= ZTOT_RCDRYG(:) +ZRCDRYG(:) *ZMAXTIME(:) + ZTOT_RIDRYG(:)= ZTOT_RIDRYG(:) +ZRIDRYG(:) *ZMAXTIME(:) + ZTOT_RRDRYG(:)= ZTOT_RRDRYG(:) +ZRRDRYG(:) *ZMAXTIME(:) + ZTOT_RSDRYG(:)= ZTOT_RSDRYG(:) +ZRSDRYG(:) *ZMAXTIME(:) + ZTOT_RGMLTR(:)= ZTOT_RGMLTR(:) +ZRGMLTR(:) *ZMAXTIME(:) + ZTOT_RCWETH(:)= ZTOT_RCWETH(:) +ZRCWETH(:) *ZMAXTIME(:) + ZTOT_RIWETH(:)= ZTOT_RIWETH(:) +ZRIWETH(:) *ZMAXTIME(:) + ZTOT_RSWETH(:)= ZTOT_RSWETH(:) +ZRSWETH(:) *ZMAXTIME(:) + ZTOT_RGWETH(:)= ZTOT_RGWETH(:) +ZRGWETH(:) *ZMAXTIME(:) + ZTOT_RRWETH(:)= ZTOT_RRWETH(:) +ZRRWETH(:) *ZMAXTIME(:) + ZTOT_RCDRYH(:)= ZTOT_RCDRYH(:) +ZRCDRYH(:) *ZMAXTIME(:) + ZTOT_RIDRYH(:)= ZTOT_RIDRYH(:) +ZRIDRYH(:) *ZMAXTIME(:) + ZTOT_RSDRYH(:)= ZTOT_RSDRYH(:) +ZRSDRYH(:) *ZMAXTIME(:) + ZTOT_RRDRYH(:)= ZTOT_RRDRYH(:) +ZRRDRYH(:) *ZMAXTIME(:) + ZTOT_RGDRYH(:)= ZTOT_RGDRYH(:) +ZRGDRYH(:) *ZMAXTIME(:) + ZTOT_RDRYHG(:)= ZTOT_RDRYHG(:) +ZRDRYHG(:) *ZMAXTIME(:) + ZTOT_RHMLTR(:)= ZTOT_RHMLTR(:) +ZRHMLTR(:) *ZMAXTIME(:) + ZTOT_RIMLTC(:)= ZTOT_RIMLTC(:) +ZRIMLTC_MR(:) + ZTOT_RCBERI(:)= ZTOT_RCBERI(:) +ZRCBERI(:) *ZMAXTIME(:) + ENDIF + ! + !*** 4.5 Next loop + ! + LSOFT=.TRUE. ! We try to adjust tendencies (inner while loop) + ZTIME(:)=ZTIME(:)+ZMAXTIME(:) + ENDDO +ENDDO + + + + + + + + + +!------------------------------------------------------------------------------- +! +!* 5. UNPACKING DIAGNOSTICS +! --------------------- +! +IF(IMICRO>0) THEN + ZW(:,:,:) = 0. + ZHLC_HCF3D(:,:,:) = UNPACK(ZHLC_HCF(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:)) + ZW(:,:,:) = 0. + ZHLC_LCF3D(:,:,:) = UNPACK(ZHLC_LCF(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:)) + ZW(:,:,:) = 0. + ZHLC_HRC3D(:,:,:) = UNPACK(ZHLC_HRC(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:)) + ZW(:,:,:) = 0. + ZHLC_LRC3D(:,:,:) = UNPACK(ZHLC_LRC(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:)) + PCIT(:,:,:) = UNPACK(ZCIT(:), MASK=LDMICRO(:,:,:), FIELD=PCIT(:,:,:)) +ELSE + ZRAINFR(:,:,:)=0. + ZHLC_HCF3D(:,:,:)=0. + ZHLC_LCF3D(:,:,:)=0. + ZHLC_HRC3D(:,:,:)=0. + ZHLC_LRC3D(:,:,:)=0. + PCIT(:,:,:) = 0. +ENDIF +IF(OWARM) THEN + ZW(:,:,:)=0. + PEVAP3D(:,:,:)=UNPACK(ZRREVAV(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:)) +ENDIF +! +! +!* 6. COMPUTES THE SLOW COLD PROCESS SOURCES OUTSIDE OF LDMICRO POINTS +! ---------------------------------------------------------------- +! +CALL ICE4_NUCLEATION_WRAPPER(IIT, IJT, IKT, .NOT. LDMICRO, & + PTHT, PPABST, PRHODREF, PEXN, ZZ_LSFACT/PEXN, ZT, & + PRVT, & + PCIT, ZZ_RVHENI_MR) +ZZ_LSFACT(:,:,:)=ZZ_LSFACT(:,:,:)/PEXNREF(:,:,:) +ZZ_LVFACT(:,:,:)=ZZ_LVFACT(:,:,:)/PEXNREF(:,:,:) +ZZ_RVHENI(:,:,:) = MIN(PRVS(:,:,:), ZZ_RVHENI_MR(:,:,:)/PTSTEP) +PRIS(:,:,:)=PRIS(:,:,:)+ZZ_RVHENI(:,:,:) +PRVS(:,:,:)=PRVS(:,:,:)-ZZ_RVHENI(:,:,:) +PTHS(:,:,:)=PTHS(:,:,:) + ZZ_RVHENI(:,:,:)*ZZ_LSFACT(:,:,:) +!------------------------------------------------------------------------------- +! +!* 7. UNPACKING AND TOTAL TENDENCIES +! ------------------------------ +! +! +!*** 7.1 total tendencies limited by available species +! +! ZW_??S variables will contain the new S variables values +! +IF(LEXT_TEND) THEN + !Z..T variables contain the exeternal tendency, we substract it + ZRVT(:) = ZRVT(:) - ZEXT_RV(:) * PTSTEP + ZRCT(:) = ZRCT(:) - ZEXT_RC(:) * PTSTEP + ZRRT(:) = ZRRT(:) - ZEXT_RR(:) * PTSTEP + ZRIT(:) = ZRIT(:) - ZEXT_RI(:) * PTSTEP + ZRST(:) = ZRST(:) - ZEXT_RS(:) * PTSTEP + ZRGT(:) = ZRGT(:) - ZEXT_RG(:) * PTSTEP + IF (KRR==7) ZRHT(:) = ZRHT(:) - ZEXT_RH(:) * PTSTEP + ZTHT(:) = ZTHT(:) - ZEXT_TH(:) * PTSTEP +ENDIF +!Tendencies computed from difference between old state and new state (can be negative) +ZW_RVS(:,:,:) = (UNPACK(ZRVT(:), MASK=LDMICRO(:,:,:), FIELD=PRVT(:,:,:)) - PRVT(:,:,:))*ZINV_TSTEP +ZW_RCS(:,:,:) = (UNPACK(ZRCT(:), MASK=LDMICRO(:,:,:), FIELD=PRCT(:,:,:)) - PRCT(:,:,:))*ZINV_TSTEP +ZW_RRS(:,:,:) = (UNPACK(ZRRT(:), MASK=LDMICRO(:,:,:), FIELD=PRRT(:,:,:)) - PRRT(:,:,:))*ZINV_TSTEP +ZW_RIS(:,:,:) = (UNPACK(ZRIT(:), MASK=LDMICRO(:,:,:), FIELD=PRIT(:,:,:)) - PRIT(:,:,:))*ZINV_TSTEP +ZW_RSS(:,:,:) = (UNPACK(ZRST(:), MASK=LDMICRO(:,:,:), FIELD=PRST(:,:,:)) - PRST(:,:,:))*ZINV_TSTEP +ZW_RGS(:,:,:) = (UNPACK(ZRGT(:), MASK=LDMICRO(:,:,:), FIELD=PRGT(:,:,:)) - PRGT(:,:,:))*ZINV_TSTEP +IF(KRR==7) THEN + ZW_RHS(:,:,:) = (UNPACK(ZRHT(:), MASK=LDMICRO(:,:,:), FIELD=PRHT(:,:,:)) - PRHT(:,:,:))*ZINV_TSTEP +ELSE + ZW_RHS(:,:,:) = 0. +ENDIF +ZW_THS(:,:,:) = (ZW_RCS(:,:,:)+ZW_RRS(:,:,:) )*ZZ_LVFACT(:,:,:) + & + & (ZW_RIS(:,:,:)+ZW_RSS(:,:,:)+ZW_RGS(:,:,:)+ZW_RHS(:,:,:))*ZZ_LSFACT(:,:,:) +!We apply these tendencies to the S variables +ZW_RVS(:,:,:) = PRVS(:,:,:) + ZW_RVS(:,:,:) +ZW_RCS(:,:,:) = PRCS(:,:,:) + ZW_RCS(:,:,:) +ZW_RRS(:,:,:) = PRRS(:,:,:) + ZW_RRS(:,:,:) +ZW_RIS(:,:,:) = PRIS(:,:,:) + ZW_RIS(:,:,:) +ZW_RSS(:,:,:) = PRSS(:,:,:) + ZW_RSS(:,:,:) +ZW_RGS(:,:,:) = PRGS(:,:,:) + ZW_RGS(:,:,:) +IF(KRR==7) ZW_RHS(:,:,:) = PRHS(:,:,:) + ZW_RHS(:,:,:) +ZW_THS(:,:,:) = PTHS(:,:,:) + ZW_THS(:,:,:) +!We correct negativities with conservation +CALL CORRECT_NEGATIVITIES(KRR, ZW_RVS, ZW_RCS, ZW_RRS, & + &ZW_RIS, ZW_RSS, ZW_RGS, & + &ZW_THS, ZZ_LVFACT, ZZ_LSFACT, ZW_RHS) +! +!*** 7.2 LBU_ENABLE case +! +IF(LBU_ENABLE) THEN + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RVHENI(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRIS(:,:,:) = PRIS(:,:,:) + ZW(:,:,:) + PRVS(:,:,:) = PRVS(:,:,:) - ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*ZZ_LSFACT(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'HENU_BU_RTH') + IF (LBUDGET_RV) CALL BUDGET(PRVS(:,:,:)*PRHODJ(:,:,:), 6, 'HENU_BU_RRV') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'HENU_BU_RRI') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCHONI(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRIS(:,:,:) = PRIS(:,:,:) + ZW(:,:,:) + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'HON_BU_RTH') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'HON_BU_RRC') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'HON_BU_RRI') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RRHONG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'SFR_BU_RTH') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'SFR_BU_RRR') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'SFR_BU_RRG') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RVDEPS(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:) + PRVS(:,:,:) = PRVS(:,:,:) - ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*ZZ_LSFACT(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'DEPS_BU_RTH') + IF (LBUDGET_RV) CALL BUDGET(PRVS(:,:,:)*PRHODJ(:,:,:), 6, 'DEPS_BU_RRV') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'DEPS_BU_RRS') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RIAGGS(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:) + PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:) + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'AGGS_BU_RRI') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'AGGS_BU_RRS') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RIAUTS(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:) + PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:) + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'AUTS_BU_RRI') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'AUTS_BU_RRS') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RVDEPG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PRVS(:,:,:) = PRVS(:,:,:) - ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*ZZ_LSFACT(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'DEPG_BU_RTH') + IF (LBUDGET_RV) CALL BUDGET(PRVS(:,:,:)*PRHODJ(:,:,:), 6, 'DEPG_BU_RRV') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'DEPG_BU_RRG') + + IF(OWARM) THEN + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCAUTR(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'AUTO_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'AUTO_BU_RRR') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCACCR(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'ACCR_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'ACCR_BU_RRR') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RREVAV(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:) + PRVS(:,:,:) = PRVS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) - ZW(:,:,:)*ZZ_LVFACT(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'REVA_BU_RTH') + IF (LBUDGET_RV) CALL BUDGET(PRVS(:,:,:)*PRHODJ(:,:,:), 6, 'REVA_BU_RRV') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'REVA_BU_RRR') + ENDIF + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCRIMSS(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCRIMSG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RSRIMCG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'RIM_BU_RTH') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'RIM_BU_RRC') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'RIM_BU_RRS') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'RIM_BU_RRG') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RRACCSS(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:) + PRSS(:,:,:) = PRSS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RRACCSG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RSACCRG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'ACC_BU_RTH') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'ACC_BU_RRR') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'ACC_BU_RRS') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'ACC_BU_RRG') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RSMLTG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCMLTSR, MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'CMEL_BU_RRS') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'CMEL_BU_RRG') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'CMEL_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'CMEL_BU_RRR') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RICFRRG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RRCFRIG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RICFRR(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:) + PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'CFRZ_BU_RTH') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'CFRZ_BU_RRR') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'CFRZ_BU_RRI') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'CFRZ_BU_RRG') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCWETG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RRWETG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RIWETG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RSWETG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'WETG_BU_RTH') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'WETG_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'WETG_BU_RRR') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'WETG_BU_RRI') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'WETG_BU_RRS') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'WETG_BU_RRG') + + IF(KRR==7) THEN + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RWETGH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'GHCV_BU_RRG') + IF (LBUDGET_RH) CALL BUDGET(PRHS(:,:,:)*PRHODJ(:,:,:), 12,'GHCV_BU_RRH') + END IF + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCDRYG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RRDRYG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RIDRYG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RSDRYG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'DRYG_BU_RTH') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'DRYG_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'DRYG_BU_RRR') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'DRYG_BU_RRI') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'DRYG_BU_RRS') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'DRYG_BU_RRG') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RGMLTR(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) - ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'GMLT_BU_RTH') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'GMLT_BU_RRR') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'GMLT_BU_RRG') + + IF(KRR==7) THEN + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCWETH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RRWETH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RIWETH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RSWETH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RGWETH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'WETH_BU_RTH') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'WETH_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'WETH_BU_RRR') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'WETH_BU_RRI') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'WETH_BU_RRS') + IF (LBUDGET_RH) CALL BUDGET(PRHS(:,:,:)*PRHODJ(:,:,:), 12,'WETH_BU_RRH') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RGWETH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'HGCV_BU_RRG') + IF (LBUDGET_RH) CALL BUDGET(PRHS(:,:,:)*PRHODJ(:,:,:), 12,'HGCV_BU_RRH') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCDRYH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RRDRYH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RIDRYH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RSDRYH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRSS(:,:,:) = PRSS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RGDRYH(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRGS(:,:,:) = PRGS(:,:,:) - ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) + ZW(:,:,:) + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RDRYHG(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRHS(:,:,:) = PRHS(:,:,:) - ZW(:,:,:) + PRGS(:,:,:) = PRGS(:,:,:) + ZW(:,:,:) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'DRYH_BU_RTH') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'DRYH_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'DRYH_BU_RRR') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'DRYH_BU_RRI') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'DRYH_BU_RRS') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'DRYH_BU_RRG') + IF (LBUDGET_RH) CALL BUDGET(PRHS(:,:,:)*PRHODJ(:,:,:), 12,'DRYH_BU_RRH') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RHMLTR(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:) + PRHS(:,:,:) = PRHS(:,:,:) - ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) - ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'HMLT_BU_RTH') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'HMLT_BU_RRR') + IF (LBUDGET_RH) CALL BUDGET(PRHS(:,:,:)*PRHODJ(:,:,:), 12,'HMLT_BU_RRH') + ENDIF + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RIMLTC(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRIS(:,:,:) = PRIS(:,:,:) - ZW(:,:,:) + PRCS(:,:,:) = PRCS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) - ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'IMLT_BU_RTH') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'IMLT_BU_RRC') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'IMLT_BU_RRI') + + ZW(:,:,:) = 0. + ZW(:,:,:)=UNPACK(ZTOT_RCBERI(:), MASK=LDMICRO(:,:,:), FIELD=ZW(:,:,:))*ZINV_TSTEP + PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:) + PRIS(:,:,:) = PRIS(:,:,:) + ZW(:,:,:) + PTHS(:,:,:) = PTHS(:,:,:) + ZW(:,:,:)*(ZZ_LSFACT(:,:,:)-ZZ_LVFACT(:,:,:)) + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'BERFI_BU_RTH') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'BERFI_BU_RRC') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'BERFI_BU_RRI') +ENDIF +! +!*** 7.3 Final tendencies +! +PRVS(:,:,:) = ZW_RVS(:,:,:) +PRCS(:,:,:) = ZW_RCS(:,:,:) +PRRS(:,:,:) = ZW_RRS(:,:,:) +PRIS(:,:,:) = ZW_RIS(:,:,:) +PRSS(:,:,:) = ZW_RSS(:,:,:) +PRGS(:,:,:) = ZW_RGS(:,:,:) +IF (KRR==7) THEN + PRHS(:,:,:) = ZW_RHS(:,:,:) +ENDIF +PTHS(:,:,:) = ZW_THS(:,:,:) +IF(LBU_ENABLE) THEN + IF (LBUDGET_TH) CALL BUDGET(PTHS(:,:,:)*PRHODJ(:,:,:), 4, 'CORR_BU_RTH') + IF (LBUDGET_RV) CALL BUDGET(PRVS(:,:,:)*PRHODJ(:,:,:), 6, 'CORR_BU_RRV') + IF (LBUDGET_RC) CALL BUDGET(PRCS(:,:,:)*PRHODJ(:,:,:), 7, 'CORR_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET(PRRS(:,:,:)*PRHODJ(:,:,:), 8, 'CORR_BU_RRR') + IF (LBUDGET_RI) CALL BUDGET(PRIS(:,:,:)*PRHODJ(:,:,:), 9, 'CORR_BU_RRI') + IF (LBUDGET_RS) CALL BUDGET(PRSS(:,:,:)*PRHODJ(:,:,:), 10,'CORR_BU_RRS') + IF (LBUDGET_RG) CALL BUDGET(PRGS(:,:,:)*PRHODJ(:,:,:), 11,'CORR_BU_RRG') + IF (KRR==7) THEN + IF (LBUDGET_RH) CALL BUDGET(PRHS(:,:,:)*PRHODJ(:,:,:), 12,'CORR_BU_RRH') + ENDIF +ENDIF +! +!------------------------------------------------------------------------------- +! +!* 8. COMPUTE THE SEDIMENTATION (RS) SOURCE +! ------------------------------------- +! +IF(LSEDIM_AFTER) THEN + ! + !* 8.1 sedimentation + ! + IF(HSEDIM=='STAT') THEN + !SR: It *seems* that we must have two separate calls for ifort + IF(KRR==7) THEN + CALL ICE4_SEDIMENTATION_STAT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKTB, IKTE, IKT, KKL, & + &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, & + &PRHODREF, PPABST, PTHT, PRHODJ, & + &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,& + &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,& + &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, & + &PSEA, PTOWN, & + &PINPRH=PINPRH, PRHT=PRHS*PTSTEP, PRHS=PRHS, PFPR=PFPR) + ELSE + CALL ICE4_SEDIMENTATION_STAT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKTB, IKTE, IKT, KKL, & + &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ,& + &PRHODREF, PPABST, PTHT, PRHODJ, & + &PRCS, PRCS*PTSTEP, PRRS, PRRS*PTSTEP, PRIS, PRIS*PTSTEP,& + &PRSS, PRSS*PTSTEP, PRGS, PRGS*PTSTEP,& + &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, & + &PSEA, PTOWN, & + &PFPR=PFPR) + ENDIF + PINPRS(:,:) = PINPRS(:,:) + ZINPRI(:,:) + !No negativity correction here as we apply sedimentation on PR.S*PTSTEP variables + ELSEIF(HSEDIM=='SPLI') THEN + !SR: It *seems* that we must have two separate calls for ifort + IF(KRR==7) THEN + CALL ICE4_SEDIMENTATION_SPLIT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKTB, IKTE, IKT, KKL, & + &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, & + &PRHODREF, PPABST, PTHT, PRHODJ, & + &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,& + &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, & + &PSEA, PTOWN, & + &PINPRH=PINPRH, PRHT=PRHT, PRHS=PRHS, PFPR=PFPR) + ELSE + CALL ICE4_SEDIMENTATION_SPLIT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKTB, IKTE, IKT, KKL, & + &PTSTEP, KRR, OSEDIC, LDEPOSC, XVDEPOSC, PDZZ, & + &PRHODREF, PPABST, PTHT, PRHODJ, & + &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,& + &PINPRC, PINDEP, PINPRR, ZINPRI, PINPRS, PINPRG, & + &PSEA, PTOWN, & + &PFPR=PFPR) + ENDIF + PINPRS(:,:) = PINPRS(:,:) + ZINPRI(:,:) + !We correct negativities with conservation + !SPLI algorith uses a time-splitting. Inside the loop a temporary m.r. is used. + ! It is initialized with the m.r. at T and is modified by two tendencies: + ! sedimentation tendency and an external tendency which represents all other + ! processes (mainly advection and microphysical processes). If both tendencies + ! are negative, sedimentation can remove a specie at a given sub-timestep. From + ! this point sedimentation stops for the remaining sub-timesteps but the other tendency + ! will be still active and will lead to negative values. + ! We could prevent the algorithm to not consume too much a specie, instead we apply + ! a correction here. + CALL CORRECT_NEGATIVITIES(KRR, PRVS, PRCS, PRRS, & + &PRIS, PRSS, PRGS, & + &PTHS, ZZ_LVFACT, ZZ_LSFACT, PRHS) + ELSE + WRITE(*,*) ' STOP' + WRITE(*,*) ' NO SEDIMENTATION SCHEME FOR HSEDIM=', HSEDIM + CALL PRINT_MSG(NVERB_FATAL,'GEN','RAIN_ICE_RED','') + END IF + ! + !* 8.2 budget storage + ! + IF (LBUDGET_RC .AND. OSEDIC) & + CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:), 7 , 'SEDI_BU_RRC') + IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:), 8 , 'SEDI_BU_RRR') + IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:), 9 , 'SEDI_BU_RRI') + IF (LBUDGET_RS) CALL BUDGET (PRSS(:,:,:)*PRHODJ(:,:,:), 10, 'SEDI_BU_RRS') + IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:), 11, 'SEDI_BU_RRG') + IF ( KRR == 7 .AND. LBUDGET_RH) & + CALL BUDGET (PRHS(:,:,:)*PRHODJ(:,:,:), 12, 'SEDI_BU_RRH') + ! + !sedimentation of rain fraction + CALL ICE4_RAINFR_VERT(IIB, IIE, IIT, IJB, IJE, IJT, IKB, IKE, IKT, KKL, ZRAINFR, PRRS(:,:,:)*PTSTEP) + +ENDIF +! +! +CONTAINS + FUNCTION RAIN_ICE_COUNTJV(LTAB, KIT, KJT, KKT, KSIZE, I1,I2,I3) RESULT(IC) + ! + !* 0. DECLARATIONS + ! ------------ + ! + IMPLICIT NONE + ! + !* 0.2 declaration of local variables + ! + ! + INTEGER, INTENT(IN) :: KIT, KJT, KKT, KSIZE + LOGICAL, DIMENSION(KIT,KJT,KKT), INTENT(IN) :: LTAB ! Mask + INTEGER, DIMENSION(KSIZE), INTENT(OUT) :: I1,I2,I3 ! Used to replace the COUNT and PACK + INTEGER :: JI,JJ,JK,IC + ! + !------------------------------------------------------------------------------- + ! + IC = 0 + DO JK = 1, SIZE(LTAB,3) + DO JJ = 1, SIZE(LTAB,2) + DO JI = 1, SIZE(LTAB,1) + IF(LTAB(JI,JJ,JK)) THEN + IC = IC +1 + I1(IC) = JI + I2(IC) = JJ + I3(IC) = JK + END IF + END DO + END DO + END DO + ! + ! + END FUNCTION RAIN_ICE_COUNTJV + ! + ! + SUBROUTINE CORRECT_NEGATIVITIES(KRR, PRV, PRC, PRR, & + &PRI, PRS, PRG, & + &PTH, PLVFACT, PLSFACT, PRH) + ! + IMPLICIT NONE + ! + INTEGER, INTENT(IN) :: KRR + REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRV, PRC, PRR, PRI, PRS, PRG, PTH + REAL, DIMENSION(:,:,:), INTENT(IN) :: PLVFACT, PLSFACT + REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: PRH + ! + REAL, DIMENSION(SIZE(PRV,1), SIZE(PRV,2), SIZE(PRV,3)) :: ZW + ! + !We correct negativities with conservation + ! 1) deal with negative values for mixing ratio, except for vapor + WHERE(PRC(:,:,:)<0.) + PRV(:,:,:)=PRV(:,:,:)+PRC(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-PRC(:,:,:)*PLVFACT(:,:,:) + PRC(:,:,:)=0. + ENDWHERE + WHERE(PRR(:,:,:)<0.) + PRV(:,:,:)=PRV(:,:,:)+PRR(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-PRR(:,:,:)*PLVFACT(:,:,:) + PRR(:,:,:)=0. + ENDWHERE + WHERE(PRI(:,:,:)<0.) + PRV(:,:,:)=PRV(:,:,:)+PRI(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-PRI(:,:,:)*PLSFACT(:,:,:) + PRI(:,:,:)=0. + ENDWHERE + WHERE(PRS(:,:,:)<0.) + PRV(:,:,:)=PRV(:,:,:)+PRS(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-PRS(:,:,:)*PLSFACT(:,:,:) + PRS(:,:,:)=0. + ENDWHERE + WHERE(PRG(:,:,:)<0.) + PRV(:,:,:)=PRV(:,:,:)+PRG(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-PRG(:,:,:)*PLSFACT(:,:,:) + PRG(:,:,:)=0. + ENDWHERE + IF(KRR==7) THEN + WHERE(PRH(:,:,:)<0.) + PRV(:,:,:)=PRV(:,:,:)+PRH(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-PRH(:,:,:)*PLSFACT(:,:,:) + PRH(:,:,:)=0. + ENDWHERE + ENDIF + ! 2) deal with negative vapor mixing ratio + WHERE(PRV(:,:,:)<0. .AND. PRC(:,:,:)+PRI(:,:,:)>0.) + ! for rc and ri, we keep ice fraction constant + ZW(:,:,:)=MIN(1., -PRV(:,:,:)/(PRC(:,:,:)+PRI(:,:,:))) ! Proportion of rc+ri to convert into rv + PTH(:,:,:)=PTH(:,:,:)-ZW(:,:,:)*(PRC(:,:,:)*PLVFACT(:,:,:)+PRI(:,:,:)*PLSFACT(:,:,:)) + PRV(:,:,:)=PRV(:,:,:)+ZW(:,:,:)*(PRC(:,:,:)+PRI(:,:,:)) + PRC(:,:,:)=(1.-ZW(:,:,:))*PRC(:,:,:) + PRI(:,:,:)=(1.-ZW(:,:,:))*PRI(:,:,:) + ENDWHERE + WHERE(PRV(:,:,:)<0. .AND. PRR(:,:,:)>0.) + ZW(:,:,:)=MIN(PRR(:,:,:), -PRV(:,:,:)) ! Quantity of rr to convert into rv + PRV(:,:,:)=PRV(:,:,:)+ZW(:,:,:) + PRR(:,:,:)=PRR(:,:,:)-ZW(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-ZW(:,:,:)*PLVFACT(:,:,:) + ENDWHERE + WHERE(PRV(:,:,:)<0. .AND. PRS(:,:,:)>0.) + ZW(:,:,:)=MIN(PRS(:,:,:), -PRV(:,:,:)) ! Quantity of rs to convert into rv + PRV(:,:,:)=PRV(:,:,:)+ZW(:,:,:) + PRS(:,:,:)=PRS(:,:,:)-ZW(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-ZW(:,:,:)*PLSFACT(:,:,:) + ENDWHERE + WHERE(PRV(:,:,:)<0. .AND. PRG(:,:,:)>0.) + ZW(:,:,:)=MIN(PRG(:,:,:), -PRV(:,:,:)) ! Quantity of rg to convert into rv + PRV(:,:,:)=PRV(:,:,:)+ZW(:,:,:) + PRG(:,:,:)=PRG(:,:,:)-ZW(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-ZW(:,:,:)*PLSFACT(:,:,:) + ENDWHERE + IF(KRR==7) THEN + WHERE(PRV(:,:,:)<0. .AND. PRH(:,:,:)>0.) + ZW(:,:,:)=MIN(PRH(:,:,:), -PRV(:,:,:)) ! Quantity of rh to convert into rv + PRV(:,:,:)=PRV(:,:,:)+ZW(:,:,:) + PRH(:,:,:)=PRH(:,:,:)-ZW(:,:,:) + PTH(:,:,:)=PTH(:,:,:)-ZW(:,:,:)*PLSFACT(:,:,:) + ENDWHERE + ENDIF + ! + ! + END SUBROUTINE CORRECT_NEGATIVITIES +! + recursive subroutine quicksort(a, first, last) + implicit none + integer a(*), x, t + integer first, last + integer i, j + + x = a( (first+last) / 2 ) + i = first + j = last + do + do while (a(i) < x) + i=i+1 + end do + do while (x < a(j)) + j=j-1 + end do + if (i >= j) exit + t = a(i); a(i) = a(j); a(j) = t + i=i+1 + j=j-1 + end do + if (first < i-1) call quicksort(a, first, i-1) + if (j+1 < last) call quicksort(a, j+1, last) + end subroutine quicksort + +END SUBROUTINE RAIN_ICE_RED -- GitLab