| Line |
Branch |
Exec |
Source |
| 1 |
|
|
! |
| 2 |
|
|
! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $ |
| 3 |
|
|
! |
| 4 |
|
✗ |
SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, & |
| 5 |
|
✗ |
& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & |
| 6 |
|
✗ |
& lalim,f0,detr_star,entr_star,f_star,csc,ztva, & |
| 7 |
|
✗ |
& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & |
| 8 |
|
|
& ,lev_out,lunout1,igout) |
| 9 |
|
|
! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam) |
| 10 |
|
|
!-------------------------------------------------------------------------- |
| 11 |
|
|
! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam |
| 12 |
|
|
! |
| 13 |
|
|
!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance |
| 14 |
|
|
! This versions starts from a cleaning of thermcell_plume_6A (2019/01/20) |
| 15 |
|
|
! thermcell_plume_6A is activate for flag_thermas_ed < 10 |
| 16 |
|
|
! thermcell_plume_5B for flag_thermas_ed < 20 |
| 17 |
|
|
! thermcell_plume for flag_thermals_ed>= 20 |
| 18 |
|
|
! Various options are controled by the flag_thermals_ed parameter |
| 19 |
|
|
! = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8 |
| 20 |
|
|
! = 21 : the Jam strato-cumulus modif is not activated in detrainment |
| 21 |
|
|
! = 29 : an other way to compute the modified buoyancy (to be tested) |
| 22 |
|
|
!-------------------------------------------------------------------------- |
| 23 |
|
|
USE IOIPSL, ONLY : getin |
| 24 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
| 25 |
|
|
|
| 26 |
|
|
USE print_control_mod, ONLY: prt_level |
| 27 |
|
|
IMPLICIT NONE |
| 28 |
|
|
|
| 29 |
|
|
! |
| 30 |
|
|
! $Header$ |
| 31 |
|
|
! |
| 32 |
|
|
! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre |
| 33 |
|
|
! veillez � n'utiliser que des ! pour les commentaires |
| 34 |
|
|
! et � bien positionner les & des lignes de continuation |
| 35 |
|
|
! (les placer en colonne 6 et en colonne 73) |
| 36 |
|
|
! |
| 37 |
|
|
! |
| 38 |
|
|
! A1.0 Fundamental constants |
| 39 |
|
|
REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO |
| 40 |
|
|
! A1.1 Astronomical constants |
| 41 |
|
|
REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA |
| 42 |
|
|
! A1.1.bis Constantes concernant l'orbite de la Terre: |
| 43 |
|
|
REAL R_ecc, R_peri, R_incl |
| 44 |
|
|
! A1.2 Geoide |
| 45 |
|
|
REAL RA,RG,R1SA |
| 46 |
|
|
! A1.3 Radiation |
| 47 |
|
|
! REAL RSIGMA,RI0 |
| 48 |
|
|
REAL RSIGMA |
| 49 |
|
|
! A1.4 Thermodynamic gas phase |
| 50 |
|
|
REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12 |
| 51 |
|
|
REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV |
| 52 |
|
|
REAL RKAPPA,RETV, eps_w |
| 53 |
|
|
! A1.5,6 Thermodynamic liquid,solid phases |
| 54 |
|
|
REAL RCW,RCS |
| 55 |
|
|
! A1.7 Thermodynamic transition of phase |
| 56 |
|
|
REAL RLVTT,RLSTT,RLMLT,RTT,RATM |
| 57 |
|
|
! A1.8 Curve of saturation |
| 58 |
|
|
REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS |
| 59 |
|
|
REAL RALPD,RBETD,RGAMD |
| 60 |
|
|
! |
| 61 |
|
|
COMMON/YOMCST/RPI ,RCLUM ,RHPLA ,RKBOL ,RNAVO & |
| 62 |
|
|
& ,RDAY ,REA ,REPSM ,RSIYEA,RSIDAY,ROMEGA & |
| 63 |
|
|
& ,R_ecc, R_peri, R_incl & |
| 64 |
|
|
& ,RA ,RG ,R1SA & |
| 65 |
|
|
& ,RSIGMA & |
| 66 |
|
|
& ,R ,RMD ,RMV ,RD ,RV ,RCPD & |
| 67 |
|
|
& ,RMO3 ,RMCO2 ,RMC ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12 & |
| 68 |
|
|
& ,RCPV ,RCVD ,RCVV ,RKAPPA,RETV, eps_w & |
| 69 |
|
|
& ,RCW ,RCS & |
| 70 |
|
|
& ,RLVTT ,RLSTT ,RLMLT ,RTT ,RATM & |
| 71 |
|
|
& ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS & |
| 72 |
|
|
& ,RALPD ,RBETD ,RGAMD |
| 73 |
|
|
! ------------------------------------------------------------------ |
| 74 |
|
|
!$OMP THREADPRIVATE(/YOMCST/) |
| 75 |
|
|
! |
| 76 |
|
|
! $Id: YOETHF.h 2799 2017-02-24 18:50:33Z jyg $ |
| 77 |
|
|
! |
| 78 |
|
|
! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre |
| 79 |
|
|
! veillez n'utiliser que des ! pour les commentaires |
| 80 |
|
|
! et bien positionner les & des lignes de continuation |
| 81 |
|
|
! (les placer en colonne 6 et en colonne 73) |
| 82 |
|
|
! |
| 83 |
|
|
!* COMMON *YOETHF* DERIVED CONSTANTS SPECIFIC TO ECMWF THERMODYNAMICS |
| 84 |
|
|
! |
| 85 |
|
|
! *R__ES* *CONSTANTS USED FOR COMPUTATION OF SATURATION |
| 86 |
|
|
! MIXING RATIO OVER LIQUID WATER(*R_LES*) OR |
| 87 |
|
|
! ICE(*R_IES*). |
| 88 |
|
|
! *RVTMP2* *RVTMP2=RCPV/RCPD-1. |
| 89 |
|
|
! *RHOH2O* *DENSITY OF LIQUID WATER. (RATM/100.) |
| 90 |
|
|
! |
| 91 |
|
|
REAL R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES |
| 92 |
|
|
REAL RVTMP2, RHOH2O |
| 93 |
|
|
REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU |
| 94 |
|
|
REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2 |
| 95 |
|
|
LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90 |
| 96 |
|
|
! If FALSE, then variables set by suphel.F90 |
| 97 |
|
|
COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES, & |
| 98 |
|
|
& RVTMP2, RHOH2O, & |
| 99 |
|
|
& R5ALVCP,R5ALSCP,RALVDCP,RALSDCP, & |
| 100 |
|
|
& RALFDCP,RTWAT,RTBER,RTBERCU, & |
| 101 |
|
|
& RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,& |
| 102 |
|
|
& RKOOP2, & |
| 103 |
|
|
& OK_BAD_ECMWF_THERMO |
| 104 |
|
|
|
| 105 |
|
|
!$OMP THREADPRIVATE(/YOETHF/) |
| 106 |
|
|
! |
| 107 |
|
|
! $Header$ |
| 108 |
|
|
! |
| 109 |
|
|
! |
| 110 |
|
|
! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre |
| 111 |
|
|
! veillez n'utiliser que des ! pour les commentaires |
| 112 |
|
|
! et bien positionner les & des lignes de continuation |
| 113 |
|
|
! (les placer en colonne 6 et en colonne 73) |
| 114 |
|
|
! |
| 115 |
|
|
! ------------------------------------------------------------------ |
| 116 |
|
|
! This COMDECK includes the Thermodynamical functions for the cy39 |
| 117 |
|
|
! ECMWF Physics package. |
| 118 |
|
|
! Consistent with YOMCST Basic physics constants, assuming the |
| 119 |
|
|
! partial pressure of water vapour is given by a first order |
| 120 |
|
|
! Taylor expansion of Qs(T) w.r.t. to Temperature, using constants |
| 121 |
|
|
! in YOETHF |
| 122 |
|
|
! ------------------------------------------------------------------ |
| 123 |
|
|
REAL PTARG, PDELARG, P5ARG, PQSARG, PCOARG |
| 124 |
|
|
REAL FOEEW, FOEDE, qsats, qsatl, dqsats, dqsatl |
| 125 |
|
|
LOGICAL thermcep |
| 126 |
|
|
PARAMETER (thermcep=.TRUE.) |
| 127 |
|
|
! |
| 128 |
|
|
FOEEW ( PTARG,PDELARG ) = EXP ( & |
| 129 |
|
|
& (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) & |
| 130 |
|
|
& / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) ) |
| 131 |
|
|
! |
| 132 |
|
|
FOEDE ( PTARG,PDELARG,P5ARG,PQSARG,PCOARG ) = PQSARG*PCOARG*P5ARG & |
| 133 |
|
|
& / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2 |
| 134 |
|
|
! |
| 135 |
|
|
qsats(ptarg) = 100.0 * 0.622 * 10.0 & |
| 136 |
|
|
& ** (2.07023 - 0.00320991 * ptarg & |
| 137 |
|
|
& - 2484.896 / ptarg + 3.56654 * LOG10(ptarg)) |
| 138 |
|
|
qsatl(ptarg) = 100.0 * 0.622 * 10.0 & |
| 139 |
|
|
& ** (23.8319 - 2948.964 / ptarg & |
| 140 |
|
|
& - 5.028 * LOG10(ptarg) & |
| 141 |
|
|
& - 29810.16 * EXP( - 0.0699382 * ptarg) & |
| 142 |
|
|
& + 25.21935 * EXP( - 2999.924 / ptarg)) |
| 143 |
|
|
! |
| 144 |
|
|
dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg & |
| 145 |
|
|
& +2484.896*LOG(10.)/ptarg**2 & |
| 146 |
|
|
& -0.00320991*LOG(10.)) |
| 147 |
|
|
dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* & |
| 148 |
|
|
& (2948.964/ptarg**2-5.028/LOG(10.)/ptarg & |
| 149 |
|
|
& +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) & |
| 150 |
|
|
& +29810.16*0.0699382*EXP(-0.0699382*ptarg)) |
| 151 |
|
|
integer :: iflag_thermals,nsplit_thermals |
| 152 |
|
|
|
| 153 |
|
|
!!! nrlmd le 10/04/2012 |
| 154 |
|
|
integer :: iflag_trig_bl,iflag_clos_bl |
| 155 |
|
|
integer :: tau_trig_shallow,tau_trig_deep |
| 156 |
|
|
real :: s_trig |
| 157 |
|
|
!!! fin nrlmd le 10/04/2012 |
| 158 |
|
|
|
| 159 |
|
|
real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30. |
| 160 |
|
|
real :: alp_bl_k |
| 161 |
|
|
real :: tau_thermals,fact_thermals_ed_dz |
| 162 |
|
|
integer,parameter :: w2di_thermals=0 |
| 163 |
|
|
integer :: isplit |
| 164 |
|
|
|
| 165 |
|
|
integer :: iflag_coupl,iflag_clos,iflag_wake |
| 166 |
|
|
integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure |
| 167 |
|
|
|
| 168 |
|
|
common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure |
| 169 |
|
|
common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz |
| 170 |
|
|
common/ctherm4/iflag_coupl,iflag_clos,iflag_wake |
| 171 |
|
|
common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux |
| 172 |
|
|
|
| 173 |
|
|
!!! nrlmd le 10/04/2012 |
| 174 |
|
|
common/ctherm6/iflag_trig_bl,iflag_clos_bl |
| 175 |
|
|
common/ctherm7/tau_trig_shallow,tau_trig_deep |
| 176 |
|
|
common/ctherm8/s_trig |
| 177 |
|
|
!!! fin nrlmd le 10/04/2012 |
| 178 |
|
|
|
| 179 |
|
|
!$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/) |
| 180 |
|
|
!$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/) |
| 181 |
|
|
|
| 182 |
|
|
INTEGER itap |
| 183 |
|
|
INTEGER lunout1,igout |
| 184 |
|
|
INTEGER ngrid,klev |
| 185 |
|
|
REAL ptimestep |
| 186 |
|
|
REAL ztv(ngrid,klev) |
| 187 |
|
|
REAL zthl(ngrid,klev) |
| 188 |
|
|
REAL po(ngrid,klev) |
| 189 |
|
|
REAL zl(ngrid,klev) |
| 190 |
|
|
REAL rhobarz(ngrid,klev) |
| 191 |
|
|
REAL zlev(ngrid,klev+1) |
| 192 |
|
|
REAL pplev(ngrid,klev+1) |
| 193 |
|
|
REAL pphi(ngrid,klev) |
| 194 |
|
|
REAL zpspsk(ngrid,klev) |
| 195 |
|
|
REAL alim_star(ngrid,klev) |
| 196 |
|
|
REAL f0(ngrid) |
| 197 |
|
|
INTEGER lalim(ngrid) |
| 198 |
|
|
integer lev_out ! niveau pour les print |
| 199 |
|
|
integer nbpb |
| 200 |
|
|
|
| 201 |
|
|
real alim_star_tot(ngrid) |
| 202 |
|
|
|
| 203 |
|
|
REAL ztva(ngrid,klev) |
| 204 |
|
|
REAL ztla(ngrid,klev) |
| 205 |
|
|
REAL zqla(ngrid,klev) |
| 206 |
|
|
REAL zqta(ngrid,klev) |
| 207 |
|
|
REAL zha(ngrid,klev) |
| 208 |
|
|
|
| 209 |
|
|
REAL detr_star(ngrid,klev) |
| 210 |
|
|
REAL coefc |
| 211 |
|
|
REAL entr_star(ngrid,klev) |
| 212 |
|
✗ |
REAL detr(ngrid,klev) |
| 213 |
|
✗ |
REAL entr(ngrid,klev) |
| 214 |
|
|
|
| 215 |
|
|
REAL csc(ngrid,klev) |
| 216 |
|
|
|
| 217 |
|
|
REAL zw2(ngrid,klev+1) |
| 218 |
|
|
REAL w_est(ngrid,klev+1) |
| 219 |
|
|
REAL f_star(ngrid,klev+1) |
| 220 |
|
✗ |
REAL wa_moy(ngrid,klev+1) |
| 221 |
|
|
|
| 222 |
|
|
REAL ztva_est(ngrid,klev) |
| 223 |
|
✗ |
REAL ztv_est(ngrid,klev) |
| 224 |
|
✗ |
REAL zqla_est(ngrid,klev) |
| 225 |
|
|
REAL zqsatth(ngrid,klev) |
| 226 |
|
✗ |
REAL zta_est(ngrid,klev) |
| 227 |
|
✗ |
REAL ztemp(ngrid),zqsat(ngrid) |
| 228 |
|
|
REAL zdw2,zdw2bis |
| 229 |
|
|
REAL zw2modif |
| 230 |
|
|
REAL zw2fact,zw2factbis |
| 231 |
|
✗ |
REAL zeps(ngrid,klev) |
| 232 |
|
|
|
| 233 |
|
|
REAL linter(ngrid) |
| 234 |
|
|
INTEGER lmix(ngrid) |
| 235 |
|
|
INTEGER lmix_bis(ngrid) |
| 236 |
|
✗ |
REAL wmaxa(ngrid) |
| 237 |
|
|
|
| 238 |
|
|
INTEGER ig,l,k,lt,it,lm |
| 239 |
|
|
|
| 240 |
|
✗ |
real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m |
| 241 |
|
✗ |
real zbuoyjam(ngrid,klev),zdqtjam(ngrid,klev) |
| 242 |
|
|
real zdz2,zdz3,lmel,entrbis,zdzbis |
| 243 |
|
✗ |
real d_temp(ngrid) |
| 244 |
|
|
real ztv1,ztv2,factinv,zinv,zlmel |
| 245 |
|
|
real zlmelup,zlmeldwn,zlt,zltdwn,zltup |
| 246 |
|
|
real atv1,atv2,btv1,btv2 |
| 247 |
|
|
real ztv_est1,ztv_est2 |
| 248 |
|
|
real zcor,zdelta,zcvm5,qlbef |
| 249 |
|
|
real zbetalpha, coefzlmel |
| 250 |
|
|
real eps |
| 251 |
|
|
REAL REPS,RLvCp,DDT0 |
| 252 |
|
|
PARAMETER (DDT0=.01) |
| 253 |
|
|
logical Zsat |
| 254 |
|
✗ |
LOGICAL active(ngrid),activetmp(ngrid) |
| 255 |
|
|
REAL fact_gamma,fact_gamma2,fact_epsilon2 |
| 256 |
|
|
|
| 257 |
|
|
REAL, SAVE :: fact_epsilon=0.002 |
| 258 |
|
|
REAL, SAVE :: betalpha=0.9 |
| 259 |
|
|
REAL, SAVE :: afact=2./3. |
| 260 |
|
|
REAL, SAVE :: fact_shell=1. |
| 261 |
|
|
REAL,SAVE :: detr_min=1.e-5 |
| 262 |
|
|
REAL,SAVE :: entr_min=1.e-5 |
| 263 |
|
|
REAL,SAVE :: detr_q_coef=0.012 |
| 264 |
|
|
REAL,SAVE :: detr_q_power=0.5 |
| 265 |
|
|
REAL,SAVE :: mix0=0. |
| 266 |
|
|
INTEGER,SAVE :: thermals_flag_alim=0 |
| 267 |
|
|
|
| 268 |
|
|
!$OMP THREADPRIVATE(fact_epsilon, betalpha, afact, fact_shell) |
| 269 |
|
|
!$OMP THREADPRIVATE(detr_min, entr_min, detr_q_coef, detr_q_power) |
| 270 |
|
|
!$OMP THREADPRIVATE( mix0, thermals_flag_alim) |
| 271 |
|
|
|
| 272 |
|
|
LOGICAL, SAVE :: first=.true. |
| 273 |
|
|
!$OMP THREADPRIVATE(first) |
| 274 |
|
|
|
| 275 |
|
|
|
| 276 |
|
|
REAL c2(ngrid,klev) |
| 277 |
|
|
|
| 278 |
|
✗ |
if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11' |
| 279 |
|
|
Zsat=.false. |
| 280 |
|
|
! Initialisation |
| 281 |
|
|
|
| 282 |
|
✗ |
RLvCp = RLVTT/RCPD |
| 283 |
|
✗ |
IF (first) THEN |
| 284 |
|
|
|
| 285 |
|
✗ |
CALL getin_p('thermals_fact_epsilon',fact_epsilon) |
| 286 |
|
✗ |
CALL getin_p('thermals_betalpha',betalpha) |
| 287 |
|
✗ |
CALL getin_p('thermals_afact',afact) |
| 288 |
|
✗ |
CALL getin_p('thermals_fact_shell',fact_shell) |
| 289 |
|
✗ |
CALL getin_p('thermals_detr_min',detr_min) |
| 290 |
|
✗ |
CALL getin_p('thermals_entr_min',entr_min) |
| 291 |
|
✗ |
CALL getin_p('thermals_detr_q_coef',detr_q_coef) |
| 292 |
|
✗ |
CALL getin_p('thermals_detr_q_power',detr_q_power) |
| 293 |
|
✗ |
CALL getin_p('thermals_mix0',mix0) |
| 294 |
|
✗ |
CALL getin_p('thermals_flag_alim',thermals_flag_alim) |
| 295 |
|
|
|
| 296 |
|
|
|
| 297 |
|
✗ |
first=.false. |
| 298 |
|
|
ENDIF |
| 299 |
|
|
|
| 300 |
|
✗ |
zbetalpha=betalpha/(1.+betalpha) |
| 301 |
|
|
|
| 302 |
|
|
|
| 303 |
|
|
! Initialisations des variables r?elles |
| 304 |
|
|
if (1==1) then |
| 305 |
|
✗ |
ztva(:,:)=ztv(:,:) |
| 306 |
|
✗ |
ztva_est(:,:)=ztva(:,:) |
| 307 |
|
✗ |
ztv_est(:,:)=ztv(:,:) |
| 308 |
|
✗ |
ztla(:,:)=zthl(:,:) |
| 309 |
|
✗ |
zqta(:,:)=po(:,:) |
| 310 |
|
✗ |
zqla(:,:)=0. |
| 311 |
|
✗ |
zha(:,:) = ztva(:,:) |
| 312 |
|
|
else |
| 313 |
|
|
ztva(:,:)=0. |
| 314 |
|
|
ztv_est(:,:)=0. |
| 315 |
|
|
ztva_est(:,:)=0. |
| 316 |
|
|
ztla(:,:)=0. |
| 317 |
|
|
zqta(:,:)=0. |
| 318 |
|
|
zha(:,:) =0. |
| 319 |
|
|
endif |
| 320 |
|
|
|
| 321 |
|
✗ |
zqla_est(:,:)=0. |
| 322 |
|
✗ |
zqsatth(:,:)=0. |
| 323 |
|
✗ |
zqla(:,:)=0. |
| 324 |
|
✗ |
detr_star(:,:)=0. |
| 325 |
|
✗ |
entr_star(:,:)=0. |
| 326 |
|
✗ |
alim_star(:,:)=0. |
| 327 |
|
✗ |
alim_star_tot(:)=0. |
| 328 |
|
✗ |
csc(:,:)=0. |
| 329 |
|
✗ |
detr(:,:)=0. |
| 330 |
|
✗ |
entr(:,:)=0. |
| 331 |
|
✗ |
zw2(:,:)=0. |
| 332 |
|
✗ |
zbuoy(:,:)=0. |
| 333 |
|
✗ |
zbuoyjam(:,:)=0. |
| 334 |
|
✗ |
gamma(:,:)=0. |
| 335 |
|
✗ |
zeps(:,:)=0. |
| 336 |
|
✗ |
w_est(:,:)=0. |
| 337 |
|
✗ |
f_star(:,:)=0. |
| 338 |
|
✗ |
wa_moy(:,:)=0. |
| 339 |
|
✗ |
linter(:)=1. |
| 340 |
|
|
! linter(:)=1. |
| 341 |
|
|
! Initialisation des variables entieres |
| 342 |
|
✗ |
lmix(:)=1 |
| 343 |
|
✗ |
lmix_bis(:)=2 |
| 344 |
|
✗ |
wmaxa(:)=0. |
| 345 |
|
|
|
| 346 |
|
|
|
| 347 |
|
|
!------------------------------------------------------------------------- |
| 348 |
|
|
! On ne considere comme actif que les colonnes dont les deux premieres |
| 349 |
|
|
! couches sont instables. |
| 350 |
|
|
!------------------------------------------------------------------------- |
| 351 |
|
|
|
| 352 |
|
✗ |
active(:)=ztv(:,1)>ztv(:,2) |
| 353 |
|
✗ |
d_temp(:)=0. ! Pour activer un contraste de temperature a la base |
| 354 |
|
|
! du panache |
| 355 |
|
|
! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main |
| 356 |
|
✗ |
CALL thermcell_alim(thermals_flag_alim,ngrid,klev,ztv,d_temp,zlev,alim_star,lalim) |
| 357 |
|
|
|
| 358 |
|
|
!------------------------------------------------------------------------------ |
| 359 |
|
|
! Calcul dans la premiere couche |
| 360 |
|
|
! On decide dans cette version que le thermique n'est actif que si la premiere |
| 361 |
|
|
! couche est instable. |
| 362 |
|
|
! Pourrait etre change si on veut que le thermiques puisse se d??clencher |
| 363 |
|
|
! dans une couche l>1 |
| 364 |
|
|
!------------------------------------------------------------------------------ |
| 365 |
|
✗ |
do ig=1,ngrid |
| 366 |
|
|
! Le panache va prendre au debut les caracteristiques de l'air contenu |
| 367 |
|
|
! dans cette couche. |
| 368 |
|
✗ |
if (active(ig)) then |
| 369 |
|
✗ |
ztla(ig,1)=zthl(ig,1) |
| 370 |
|
✗ |
zqta(ig,1)=po(ig,1) |
| 371 |
|
✗ |
zqla(ig,1)=zl(ig,1) |
| 372 |
|
|
!cr: attention, prise en compte de f*(1)=1 |
| 373 |
|
✗ |
f_star(ig,2)=alim_star(ig,1) |
| 374 |
|
|
zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & |
| 375 |
|
|
& *(zlev(ig,2)-zlev(ig,1)) & |
| 376 |
|
✗ |
& *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) |
| 377 |
|
✗ |
w_est(ig,2)=zw2(ig,2) |
| 378 |
|
|
endif |
| 379 |
|
|
enddo |
| 380 |
|
|
! |
| 381 |
|
|
|
| 382 |
|
|
!============================================================================== |
| 383 |
|
|
!boucle de calcul de la vitesse verticale dans le thermique |
| 384 |
|
|
!============================================================================== |
| 385 |
|
✗ |
do l=2,klev-1 |
| 386 |
|
|
!============================================================================== |
| 387 |
|
|
|
| 388 |
|
|
|
| 389 |
|
|
! On decide si le thermique est encore actif ou non |
| 390 |
|
|
! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test |
| 391 |
|
✗ |
do ig=1,ngrid |
| 392 |
|
|
active(ig)=active(ig) & |
| 393 |
|
|
& .and. zw2(ig,l)>1.e-10 & |
| 394 |
|
✗ |
& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 |
| 395 |
|
|
enddo |
| 396 |
|
|
|
| 397 |
|
|
|
| 398 |
|
|
|
| 399 |
|
|
!--------------------------------------------------------------------------- |
| 400 |
|
|
! calcul des proprietes thermodynamiques et de la vitesse de la couche l |
| 401 |
|
|
! sans tenir compte du detrainement et de l'entrainement dans cette |
| 402 |
|
|
! couche |
| 403 |
|
|
! C'est a dire qu'on suppose |
| 404 |
|
|
! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1) |
| 405 |
|
|
! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer |
| 406 |
|
|
! avant) a l'alimentation pour avoir un calcul plus propre |
| 407 |
|
|
!--------------------------------------------------------------------------- |
| 408 |
|
|
|
| 409 |
|
✗ |
ztemp(:)=zpspsk(:,l)*ztla(:,l-1) |
| 410 |
|
✗ |
call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:)) |
| 411 |
|
✗ |
do ig=1,ngrid |
| 412 |
|
|
! print*,'active',active(ig),ig,l |
| 413 |
|
✗ |
if(active(ig)) then |
| 414 |
|
✗ |
zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig)) |
| 415 |
|
✗ |
ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l) |
| 416 |
|
✗ |
zta_est(ig,l)=ztva_est(ig,l) |
| 417 |
|
✗ |
ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l) |
| 418 |
|
|
ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) & |
| 419 |
|
✗ |
& -zqla_est(ig,l))-zqla_est(ig,l)) |
| 420 |
|
|
|
| 421 |
|
|
|
| 422 |
|
|
!Modif AJAM |
| 423 |
|
|
|
| 424 |
|
✗ |
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
| 425 |
|
✗ |
zdz=zlev(ig,l+1)-zlev(ig,l) |
| 426 |
|
✗ |
lmel=fact_thermals_ed_dz*zlev(ig,l) |
| 427 |
|
|
! lmel=0.09*zlev(ig,l) |
| 428 |
|
✗ |
zlmel=zlev(ig,l)+lmel |
| 429 |
|
✗ |
zlmelup=zlmel+(zdz/2) |
| 430 |
|
|
zlmeldwn=zlmel-(zdz/2) |
| 431 |
|
|
|
| 432 |
|
|
lt=l+1 |
| 433 |
|
|
zlt=zlev(ig,lt) |
| 434 |
|
|
zdz3=zlev(ig,lt+1)-zlt |
| 435 |
|
|
zltdwn=zlt-zdz3/2 |
| 436 |
|
|
zltup=zlt+zdz3/2 |
| 437 |
|
|
|
| 438 |
|
|
!========================================================================= |
| 439 |
|
|
! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus |
| 440 |
|
|
!========================================================================= |
| 441 |
|
|
|
| 442 |
|
|
!-------------------------------------------------- |
| 443 |
|
|
lt=l+1 |
| 444 |
|
|
zlt=zlev(ig,lt) |
| 445 |
|
|
zdz2=zlev(ig,lt)-zlev(ig,l) |
| 446 |
|
|
|
| 447 |
|
✗ |
do while (lmel.gt.zdz2) |
| 448 |
|
✗ |
lt=lt+1 |
| 449 |
|
✗ |
zlt=zlev(ig,lt) |
| 450 |
|
✗ |
zdz2=zlt-zlev(ig,l) |
| 451 |
|
|
enddo |
| 452 |
|
✗ |
zdz3=zlev(ig,lt+1)-zlt |
| 453 |
|
✗ |
zltdwn=zlev(ig,lt)-zdz3/2 |
| 454 |
|
|
zlmelup=zlmel+(zdz/2) |
| 455 |
|
✗ |
coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz) |
| 456 |
|
|
zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- & |
| 457 |
|
|
& ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- & |
| 458 |
|
✗ |
& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l) |
| 459 |
|
|
|
| 460 |
|
|
!------------------------------------------------ |
| 461 |
|
|
!AJAM:nouveau calcul de w? |
| 462 |
|
|
!------------------------------------------------ |
| 463 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
| 464 |
|
|
zdzbis=zlev(ig,l)-zlev(ig,l-1) |
| 465 |
|
|
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
| 466 |
|
✗ |
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
| 467 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
| 468 |
|
✗ |
zdw2=afact*zbuoy(ig,l)/fact_epsilon |
| 469 |
|
|
zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon |
| 470 |
|
|
lm=Max(1,l-2) |
| 471 |
|
✗ |
w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) |
| 472 |
|
|
endif |
| 473 |
|
|
enddo |
| 474 |
|
|
|
| 475 |
|
|
|
| 476 |
|
|
!------------------------------------------------- |
| 477 |
|
|
!calcul des taux d'entrainement et de detrainement |
| 478 |
|
|
!------------------------------------------------- |
| 479 |
|
|
|
| 480 |
|
✗ |
do ig=1,ngrid |
| 481 |
|
✗ |
if (active(ig)) then |
| 482 |
|
|
|
| 483 |
|
|
! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1) |
| 484 |
|
✗ |
zw2m=w_est(ig,l+1) |
| 485 |
|
✗ |
zdz=zlev(ig,l+1)-zlev(ig,l) |
| 486 |
|
✗ |
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
| 487 |
|
✗ |
zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) |
| 488 |
|
✗ |
zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l) |
| 489 |
|
|
|
| 490 |
|
|
!========================================================================= |
| 491 |
|
|
! 4. Calcul de l'entrainement et du detrainement |
| 492 |
|
|
!========================================================================= |
| 493 |
|
|
|
| 494 |
|
|
detr_star(ig,l)=f_star(ig,l)*zdz & |
| 495 |
|
|
& *( mix0 * 0.1 / (zalpha+0.001) & |
| 496 |
|
|
& + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m & |
| 497 |
|
✗ |
& + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) |
| 498 |
|
|
|
| 499 |
|
✗ |
if ( iflag_thermals_ed == 20 ) then |
| 500 |
|
|
entr_star(ig,l)=f_star(ig,l)*zdz* ( & |
| 501 |
|
|
& mix0 * 0.1 / (zalpha+0.001) & |
| 502 |
|
|
& + zbetalpha*MAX(entr_min, & |
| 503 |
|
✗ |
& afact*zbuoyjam(ig,l)/zw2m - fact_epsilon)) |
| 504 |
|
|
else |
| 505 |
|
|
entr_star(ig,l)=f_star(ig,l)*zdz* ( & |
| 506 |
|
|
& mix0 * 0.1 / (zalpha+0.001) & |
| 507 |
|
|
& + zbetalpha*MAX(entr_min, & |
| 508 |
|
✗ |
& afact*zbuoy(ig,l)/zw2m - fact_epsilon)) |
| 509 |
|
|
endif |
| 510 |
|
|
|
| 511 |
|
|
! En dessous de lalim, on prend le max de alim_star et entr_star pour |
| 512 |
|
|
! alim_star et 0 sinon |
| 513 |
|
✗ |
if (l.lt.lalim(ig)) then |
| 514 |
|
✗ |
alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) |
| 515 |
|
✗ |
entr_star(ig,l)=0. |
| 516 |
|
|
endif |
| 517 |
|
|
f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & |
| 518 |
|
✗ |
& -detr_star(ig,l) |
| 519 |
|
|
|
| 520 |
|
|
endif |
| 521 |
|
|
enddo |
| 522 |
|
|
|
| 523 |
|
|
|
| 524 |
|
|
!============================================================================ |
| 525 |
|
|
! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique |
| 526 |
|
|
!=========================================================================== |
| 527 |
|
|
|
| 528 |
|
✗ |
activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 |
| 529 |
|
✗ |
do ig=1,ngrid |
| 530 |
|
✗ |
if (activetmp(ig)) then |
| 531 |
|
|
Zsat=.false. |
| 532 |
|
|
ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & |
| 533 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & |
| 534 |
|
✗ |
& /(f_star(ig,l+1)+detr_star(ig,l)) |
| 535 |
|
|
zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & |
| 536 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & |
| 537 |
|
✗ |
& /(f_star(ig,l+1)+detr_star(ig,l)) |
| 538 |
|
|
|
| 539 |
|
|
endif |
| 540 |
|
|
enddo |
| 541 |
|
|
|
| 542 |
|
✗ |
ztemp(:)=zpspsk(:,l)*ztla(:,l) |
| 543 |
|
✗ |
call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l)) |
| 544 |
|
✗ |
do ig=1,ngrid |
| 545 |
|
✗ |
if (activetmp(ig)) then |
| 546 |
|
|
! on ecrit de maniere conservative (sat ou non) |
| 547 |
|
|
! T = Tl +Lv/Cp ql |
| 548 |
|
✗ |
zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l)) |
| 549 |
|
✗ |
ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l) |
| 550 |
|
✗ |
ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l) |
| 551 |
|
|
!on rajoute le calcul de zha pour diagnostiques (temp potentielle) |
| 552 |
|
✗ |
zha(ig,l) = ztva(ig,l) |
| 553 |
|
|
ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & |
| 554 |
|
✗ |
& -zqla(ig,l))-zqla(ig,l)) |
| 555 |
|
✗ |
zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) |
| 556 |
|
✗ |
zdz=zlev(ig,l+1)-zlev(ig,l) |
| 557 |
|
|
zdzbis=zlev(ig,l)-zlev(ig,l-1) |
| 558 |
|
✗ |
zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) |
| 559 |
|
✗ |
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
| 560 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
| 561 |
|
✗ |
zdw2= afact*zbuoy(ig,l)/(fact_epsilon) |
| 562 |
|
|
zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon) |
| 563 |
|
✗ |
zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) |
| 564 |
|
|
endif |
| 565 |
|
|
enddo |
| 566 |
|
|
|
| 567 |
|
✗ |
if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l |
| 568 |
|
|
! |
| 569 |
|
|
!=========================================================================== |
| 570 |
|
|
! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max |
| 571 |
|
|
!=========================================================================== |
| 572 |
|
|
|
| 573 |
|
✗ |
nbpb=0 |
| 574 |
|
✗ |
do ig=1,ngrid |
| 575 |
|
✗ |
if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then |
| 576 |
|
|
! stop'On tombe sur le cas particulier de thermcell_dry' |
| 577 |
|
|
! print*,'On tombe sur le cas particulier de thermcell_plume' |
| 578 |
|
✗ |
nbpb=nbpb+1 |
| 579 |
|
✗ |
zw2(ig,l+1)=0. |
| 580 |
|
✗ |
linter(ig)=l+1 |
| 581 |
|
|
endif |
| 582 |
|
|
|
| 583 |
|
✗ |
if (zw2(ig,l+1).lt.0.) then |
| 584 |
|
|
linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & |
| 585 |
|
✗ |
& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) |
| 586 |
|
✗ |
zw2(ig,l+1)=0. |
| 587 |
|
|
!+CR:04/05/12:correction calcul linter pour calcul de zmax continu |
| 588 |
|
✗ |
elseif (f_star(ig,l+1).lt.0.) then |
| 589 |
|
|
linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) & |
| 590 |
|
✗ |
& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l)) |
| 591 |
|
✗ |
zw2(ig,l+1)=0. |
| 592 |
|
|
!fin CR:04/05/12 |
| 593 |
|
|
endif |
| 594 |
|
|
|
| 595 |
|
✗ |
wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) |
| 596 |
|
|
|
| 597 |
|
✗ |
if (wa_moy(ig,l+1).gt.wmaxa(ig)) then |
| 598 |
|
|
! lmix est le niveau de la couche ou w (wa_moy) est maximum |
| 599 |
|
|
!on rajoute le calcul de lmix_bis |
| 600 |
|
✗ |
if (zqla(ig,l).lt.1.e-10) then |
| 601 |
|
✗ |
lmix_bis(ig)=l+1 |
| 602 |
|
|
endif |
| 603 |
|
✗ |
lmix(ig)=l+1 |
| 604 |
|
✗ |
wmaxa(ig)=wa_moy(ig,l+1) |
| 605 |
|
|
endif |
| 606 |
|
|
enddo |
| 607 |
|
|
|
| 608 |
|
✗ |
if (nbpb>0) then |
| 609 |
|
✗ |
print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume' |
| 610 |
|
|
endif |
| 611 |
|
|
|
| 612 |
|
|
!========================================================================= |
| 613 |
|
|
! FIN DE LA BOUCLE VERTICALE |
| 614 |
|
|
enddo |
| 615 |
|
|
!========================================================================= |
| 616 |
|
|
|
| 617 |
|
|
!on recalcule alim_star_tot |
| 618 |
|
✗ |
do ig=1,ngrid |
| 619 |
|
✗ |
alim_star_tot(ig)=0. |
| 620 |
|
|
enddo |
| 621 |
|
✗ |
do ig=1,ngrid |
| 622 |
|
✗ |
do l=1,lalim(ig)-1 |
| 623 |
|
✗ |
alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) |
| 624 |
|
|
enddo |
| 625 |
|
|
enddo |
| 626 |
|
|
|
| 627 |
|
|
|
| 628 |
|
✗ |
if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l |
| 629 |
|
|
|
| 630 |
|
|
|
| 631 |
|
|
|
| 632 |
|
✗ |
return |
| 633 |
|
|
end |
| 634 |
|
|
|