Changeset 5255 for LMDZ6/trunk/libf/phylmdiso
- Timestamp:
- Oct 22, 2024, 3:35:01 PM (6 weeks ago)
- Location:
- LMDZ6/trunk/libf/phylmdiso
- Files:
-
- 4 edited
-
calwake.F90 (modified) (13 diffs)
-
isotopes_routines_mod.F90 (modified) (1 diff)
-
lmdz_wake.F90 (modified) (95 diffs)
-
physiq_mod.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmdiso/calwake.F90
r4783 r5255 6 6 dt_dwn, dq_dwn, m_dwn, m_up, dt_a, dq_a, wgen, & 7 7 sigd, Cin, & 8 wake_deltat, wake_deltaq, wake_s, awake_ dens,wake_dens, &8 wake_deltat, wake_deltaq, wake_s, awake_s, wake_dens, awake_dens, & 9 9 wake_dth, wake_h, & 10 10 wake_pe, wake_fip, wake_gfl, & … … 14 14 wake_omg, wake_dp_deltomg, & 15 15 wake_spread, wake_cstar, wake_d_deltat_gw, & 16 wake_ddeltat, wake_ddeltaq, wake_ds, awake_d dens,wake_ddens &16 wake_ddeltat, wake_ddeltaq, wake_ds, awake_ds, wake_ddens, awake_ddens & 17 17 #ifdef ISO 18 18 & ,xt,dxt_dwn,dxt_a & … … 67 67 ! ------------ 68 68 REAL, DIMENSION(klon, klev), INTENT (INOUT) :: wake_deltat, wake_deltaq 69 REAL, DIMENSION(klon), INTENT (INOUT) :: wake_s 70 REAL, DIMENSION(klon), INTENT (INOUT) :: awake_dens,wake_dens69 REAL, DIMENSION(klon), INTENT (INOUT) :: wake_s, awake_s 70 REAL, DIMENSION(klon), INTENT (INOUT) :: wake_dens, awake_dens 71 71 #ifdef ISO 72 72 REAL, DIMENSION(ntraciso,klon, klev), INTENT (INOUT) :: wake_deltaxt … … 88 88 REAL, DIMENSION(klon), INTENT (OUT) :: wake_cstar 89 89 REAL, DIMENSION(klon, klev), INTENT (OUT) :: wake_ddeltat, wake_ddeltaq 90 REAL, DIMENSION(klon), INTENT (OUT) :: wake_ds, awake_d dens,wake_ddens90 REAL, DIMENSION(klon), INTENT (OUT) :: wake_ds, awake_ds, wake_ddens, awake_ddens 91 91 #ifdef ISO 92 92 REAL, DIMENSION(ntraciso,klon, klev), INTENT (OUT) :: dxt_wake … … 115 115 REAL, DIMENSION(klon, klev) :: tx, qx 116 116 REAL, DIMENSION(klon) :: hw, wape, fip, gfl 117 REAL, DIMENSION(klon) :: sigmaw, a wdens,wdens117 REAL, DIMENSION(klon) :: sigmaw, asigmaw, wdens, awdens 118 118 REAL, DIMENSION(klon, klev) :: omgbdth 119 119 REAL, DIMENSION(klon, klev) :: dp_omgb … … 126 126 REAL, DIMENSION(klon, klev) :: d_deltat_gw 127 127 REAL, DIMENSION(klon, klev) :: d_deltatw, d_deltaqw 128 REAL, DIMENSION(klon) :: d_sigmaw, d_a wdens, d_wdens128 REAL, DIMENSION(klon) :: d_sigmaw, d_asigmaw, d_wdens, d_awdens 129 129 #ifdef ISO 130 130 REAL, DIMENSION(ntraciso,klon, klev) :: xte … … 142 142 143 143 IF (prt_level >= 10) THEN 144 print *, '-> calwake, wake_s, wgen input ',wake_s(1), wgen(1)144 print *, '-> calwake, wake_s, awake_s, wgen input ', wake_s(1), awake_s(1), wgen(1) 145 145 ENDIF 146 146 … … 191 191 d_deltaqw(:,:) = 0. 192 192 d_sigmaw(:) = 0. 193 d_asigmaw(:) = 0. 194 d_wdens(:) = 0. 193 195 d_awdens(:) = 0. 194 d_wdens(:) = 0.195 196 #ifdef ISO 196 197 dxtls(:,:,:) = 0. … … 227 228 228 229 DO i = 1, klon 229 sigmaw(i) = wake_s(i) 230 END DO 231 232 DO i = 1, klon 230 sigmaw(i) = wake_s(i) 231 asigmaw(i) = awake_s(i) 232 END DO 233 234 DO i = 1, klon 235 wdens(i) = max(0., wake_dens(i)) 233 236 awdens(i) = max(0., awake_dens(i)) 234 wdens(i) = max(0., wake_dens(i))235 237 END DO 236 238 … … 272 274 #endif 273 275 274 CALL wake(znatsurf, p, ph, pi, dtime, & 276 277 CALL wake(klon,klev,znatsurf, p, ph, pi, dtime, & 275 278 te, qe, omgbe, & 276 279 dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, & 277 280 sigd0, Cin, & 278 dtw, dqw, sigmaw, a wdens, wdens, &! state variables281 dtw, dqw, sigmaw, asigmaw, wdens, awdens, & ! state variables 279 282 dth, hw, wape, fip, gfl, & 280 283 dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, & 281 284 dtke, dqke, omg, dp_deltomg, spread, cstar, & 282 285 d_deltat_gw, & 283 d_deltatw, d_deltaqw, d_sigmaw, d_a wdens, d_wdens &! tendencies286 d_deltatw, d_deltaqw, d_sigmaw, d_asigmaw, d_wdens, d_awdens & ! tendencies 284 287 #ifdef ISO 285 288 , xte,dxtdwn,dxta,dxtw & … … 386 389 IF (ktopw(i)>0) THEN 387 390 wake_ds(i) = d_sigmaw(i)*dtime 391 awake_ds(i) = d_asigmaw(i)*dtime 388 392 awake_ddens(i) = d_awdens(i)*dtime 389 393 wake_ddens(i) = d_wdens(i)*dtime 390 394 ELSE 391 wake_ds(i) = -wake_s(i) 392 wake_ddens(i)= -wake_dens(i) 395 wake_ds(i) = -wake_s(i) 396 awake_ds(i) = -awake_s(i) 397 wake_ddens(i) = -wake_dens(i) 398 awake_ddens(i)= -awake_dens(i) 393 399 END IF 394 400 END DO … … 421 427 DO i = 1, klon 422 428 wake_s(i) = sigmaw(i) 429 awake_s(i) = asigmaw(i) 423 430 awake_dens(i) = awdens(i) 424 431 wake_dens(i) = wdens(i) … … 435 442 ENDIF ! (first) 436 443 !>jyg 444 IF (prt_level >= 10) THEN 445 print *, 'calwake ->, wake_s, awake_s ', wake_s(1), awake_s(1) 446 ENDIF 437 447 438 448 RETURN 439 449 END SUBROUTINE calwake 440 450 451 -
LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90
r5199 r5255 18822 18822 if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then 18823 18823 if (q.gt.ridicule) then 18824 write(*,*) 'xt,q=',xt,q18825 write(*,*) 'alpha=',alpha18826 write(*,*) 'toce,kcin,h0=',toce,kcin,h018827 write(*,*) 'RMerlivat=',RMerlivat18824 !write(*,*) 'xt,q=',xt,q 18825 !write(*,*) 'alpha=',alpha 18826 !write(*,*) 'toce,kcin,h0=',toce,kcin,h0 18827 !write(*,*) 'RMerlivat=',RMerlivat 18828 18828 call iso_verif_aberrant_encadre( xt/q, 'isotopes_routines_mod 18930b: iso_init_ideal') 18829 18829 endif -
LMDZ6/trunk/libf/phylmdiso/lmdz_wake.F90
r4594 r5255 1 1 MODULE lmdz_wake 2 ! $Id: wake.F90 3648 2020-03-16 15:36:59Z jghattas $ 2 3 ! $Id: lmdz_wake.F90 4908 2024-04-15 17:30:55Z jyg $ 4 3 5 CONTAINS 4 SUBROUTINE wake(znatsurf, p, ph, pi, dtime, & 5 te0, qe0, omgb, & 6 7 SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, & 8 tb0, qb0, omgb, & 6 9 dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, & 7 10 sigd_con, Cin, & 8 deltatw, deltaqw, sigmaw, a wdens,wdens, & ! state variables11 deltatw, deltaqw, sigmaw, asigmaw, wdens, awdens, & ! state variables 9 12 dth, hw, wape, fip, gfl, & 10 dtls, dqls, ktopw, omgbdth, dp_omgb, t u, qu, &11 dtke, dqke, omg, dp_deltomg, spread, cstar, &12 d_deltat_gw, & 13 d_deltatw2, d_deltaqw2, d_sigmaw2, d_a wdens2, d_wdens2 & ! tendencies13 dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, & 14 dtke, dqke, omg, dp_deltomg, wkspread, cstar, & 15 d_deltat_gw, & ! tendencies 16 d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2 & ! tendencies 14 17 15 18 #ifdef ISO 16 ,xt e0,dxtdwn,dxta,deltaxtw &17 ,dxtls,xt u,dxtke,d_deltaxtw2 &19 ,xtb0,dxtdwn,dxta,deltaxtw & 20 ,dxtls,xtx,dxtke,d_deltaxtw2 & 18 21 #endif 19 22 ) … … 29 32 ! ************************************************************** 30 33 31 USE ioipsl_getin_p_mod, ONLY : getin_p 32 USE dimphy 33 use mod_phys_lmdz_para 34 USE print_control_mod, ONLY: prt_level 34 35 USE lmdz_wake_ini , ONLY : wake_ini 36 USE lmdz_wake_ini , ONLY : prt_level,epsim1,RG,RD 37 USE lmdz_wake_ini , ONLY : stark, wdens_ref, coefgw, alpk, wk_pupper 38 USE lmdz_wake_ini , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl 39 USE lmdz_wake_ini , ONLY : ok_bug_gfl 40 USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensinit, wdensthreshold 41 USE lmdz_wake_ini , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc 42 USE lmdz_wake_ini , ONLY : iflag_wk_profile 43 USE lmdz_wake_ini , ONLY : smallestreal,wk_nsub 44 45 35 46 #ifdef ISO 36 47 USE infotrac_phy, ONLY : ntraciso=>ntiso … … 53 64 ! deltaqw : specific humidity difference between wake and off-wake regions 54 65 ! sigmaw : fractional area covered by wakes. 66 ! asigmaw : fractional area covered by active wakes. 55 67 ! wdens : number of wakes per unit area 68 ! awdens : number of active wakes per unit area 56 69 57 70 ! Variable de sortie : … … 71 84 ! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) 72 85 ! dp_deltomg : vertical gradient of omg (s-1) 73 ! spread : spreading term in d_t_wake and d_q_wake86 ! wkspread : spreading term in d_t_wake and d_q_wake 74 87 ! deltatw : updated temperature difference (T_w-T_u). 75 88 ! deltaqw : updated humidity difference (q_w-q_u). 76 89 ! sigmaw : updated wake fractional area. 90 ! asigmaw : updated active wake fractional area. 77 91 ! d_deltat_gw : delta T tendency due to GW 78 92 … … 80 94 81 95 ! aire : aire de la maille 82 ! t e0 : temperature dans l'environnement(K)83 ! q e0 : humidite dans l'environnement(kg/kg)96 ! tb0 : horizontal average of temperature (K) 97 ! qb0 : horizontal average of humidity (kg/kg) 84 98 ! omgb : vitesse verticale moyenne sur la maille (Pa/s) 85 99 ! dtdwn: source de chaleur due aux descentes (K/s) … … 101 115 ! Variables internes : 102 116 103 ! rhow : masse volumique de la poche froide 104 ! rho : environment density at P levels 105 ! rhoh : environment density at Ph levels 106 ! te : environment temperature | may change within 107 ! qe : environment humidity | sub-time-stepping 108 ! the : environment potential temperature 109 ! thu : potential temperature in undisturbed area 110 ! tu : temperature in undisturbed area 111 ! qu : humidity in undisturbed area 117 ! rho : mean density at P levels 118 ! rhoh : mean density at Ph levels 119 ! tb : mean temperature | may change within 120 ! qb : mean humidity | sub-time-stepping 121 ! thb : mean potential temperature 122 ! thx : potential temperature in (x) area 123 ! tx : temperature in (x) area 124 ! qx : humidity in (x) area 112 125 ! dp_omgb: vertical gradient og LS omega 113 126 ! omgbw : wake average vertical omega … … 115 128 ! omgbdq : flux of Delta_q transported by LS omega 116 129 ! dth : potential temperature diff. wake-undist. 117 ! th1 : first pot. temp. for vertical advection (=th u)130 ! th1 : first pot. temp. for vertical advection (=thx) 118 131 ! th2 : second pot. temp. for vertical advection (=thw) 119 132 ! q1 : first humidity for vertical advection 120 133 ! q2 : second humidity for vertical advection 121 ! d_deltatw : terme de redistribution pour deltatw122 ! d_deltaqw : terme de redistribution pour deltaqw123 ! deltatw0 : deltatw initial124 ! deltaqw0 : deltaqw initial134 ! d_deltatw : redistribution term for deltatw 135 ! d_deltaqw : redistribution term for deltaqw 136 ! deltatw0 : initial deltatw 137 ! deltaqw0 : initial deltaqw 125 138 ! hw0 : wake top hight (defined as the altitude at which deltatw=0) 126 139 ! amflux : horizontal mass flux through wake boundary 127 140 ! wdens_ref: initial number of wakes per unit area (3D) or per 128 ! unit length (2D), at the beginning of each time step129 ! Tgw : 1 sur la p ériode de onde de gravité130 ! Cgw : vitesse de propagation de onde de gravit é131 ! LL : distance entre 2 poches141 ! unit length (2D), at the beginning of each time step 142 ! Tgw : 1 sur la periode de onde de gravite 143 ! Cgw : vitesse de propagation de onde de gravite 144 ! LL : distance between 2 wakes 132 145 133 146 ! ------------------------------------------------------------------------- 134 ! D éclaration de variables147 ! Declaration de variables 135 148 ! ------------------------------------------------------------------------- 136 149 137 include "YOMCST.h"138 include "cvthermo.h"139 150 140 151 ! Arguments en entree 141 152 ! -------------------- 142 153 154 INTEGER, INTENT(IN) :: klon,klev 143 155 INTEGER, DIMENSION (klon), INTENT(IN) :: znatsurf 144 156 REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi … … 146 158 REAL, DIMENSION (klon, klev), INTENT(IN) :: omgb 147 159 REAL, INTENT(IN) :: dtime 148 REAL, DIMENSION (klon, klev), INTENT(IN) :: t e0, qe0160 REAL, DIMENSION (klon, klev), INTENT(IN) :: tb0, qb0 149 161 REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn 150 162 REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup … … 154 166 REAL, DIMENSION (klon), INTENT(IN) :: Cin 155 167 #ifdef ISO 156 REAL, DIMENSION (ntraciso,klon, klev), INTENT(IN) :: xt e0168 REAL, DIMENSION (ntraciso,klon, klev), INTENT(IN) :: xtb0 157 169 REAL, DIMENSION (ntraciso,klon, klev), INTENT(IN) :: dxtdwn 158 170 REAL, DIMENSION (ntraciso,klon, klev), INTENT(IN) :: dxta … … 164 176 REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw 165 177 REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw 178 REAL, DIMENSION (klon), INTENT(INOUT) :: asigmaw 179 REAL, DIMENSION (klon), INTENT(INOUT) :: wdens 166 180 REAL, DIMENSION (klon), INTENT(INOUT) :: awdens 167 REAL, DIMENSION (klon), INTENT(INOUT) :: wdens168 181 #ifdef ISO 169 182 REAL, DIMENSION (ntraciso,klon, klev), INTENT(INOUT) :: deltaxtw … … 174 187 175 188 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth 176 REAL, DIMENSION (klon, klev), INTENT(OUT) :: t u, qu189 REAL, DIMENSION (klon, klev), INTENT(OUT) :: tx, qx 177 190 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls 178 191 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke 179 REAL, DIMENSION (klon, klev), INTENT(OUT) :: spread ! unused (jyg)192 REAL, DIMENSION (klon, klev), INTENT(OUT) :: wkspread ! unused (jyg) 180 193 REAL, DIMENSION (klon, klev), INTENT(OUT) :: omgbdth, omg 181 194 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg 182 REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw183 195 REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar 184 196 INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw 185 ! Tendencies of state variables 197 ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields 198 ! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens) 199 REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw 186 200 REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2 187 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_awdens2, d_wdens2 201 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2 202 188 203 #ifdef ISO 189 REAL, DIMENSION (ntraciso,klon, klev), INTENT(OUT) :: xt u204 REAL, DIMENSION (ntraciso,klon, klev), INTENT(OUT) :: xtx 190 205 REAL, DIMENSION (ntraciso,klon, klev), INTENT(OUT) :: dxtls 191 206 REAL, DIMENSION (ntraciso,klon, klev), INTENT(OUT) :: dxtke … … 196 211 ! ------------------- 197 212 198 ! Variables à fixer 199 INTEGER, SAVE :: igout 200 !$OMP THREADPRIVATE(igout) 201 LOGICAL, SAVE :: first = .TRUE. 202 !$OMP THREADPRIVATE(first) 203 !jyg< 204 !! REAL, SAVE :: stark, wdens_ref, coefgw, alpk 205 REAL, SAVE, DIMENSION(2) :: wdens_ref 206 REAL, SAVE :: stark, coefgw, alpk 207 !>jyg 208 REAL, SAVE :: crep_upper, crep_sol 209 !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol) 210 211 REAL, SAVE :: tau_cv 212 !$OMP THREADPRIVATE(tau_cv) 213 REAL, SAVE :: rzero, aa0 ! minimal wake radius and area 214 !$OMP THREADPRIVATE(rzero, aa0) 215 216 LOGICAL, SAVE :: flag_wk_check_trgl 217 !$OMP THREADPRIVATE(flag_wk_check_trgl) 218 INTEGER, SAVE :: iflag_wk_check_trgl 219 !$OMP THREADPRIVATE(iflag_wk_check_trgl) 220 INTEGER, SAVE :: iflag_wk_pop_dyn 221 !$OMP THREADPRIVATE(iflag_wk_pop_dyn) 213 ! Variables a fixer 222 214 223 215 REAL :: delta_t_min 224 INTEGER :: nsub225 216 REAL :: dtimesub 226 REAL, SAVE :: wdensmin227 !$OMP THREADPRIVATE(wdensmin)228 REAL, SAVE :: sigmad, hwmin, wapecut, cstart229 !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)230 REAL, SAVE :: sigmaw_max231 !$OMP THREADPRIVATE(sigmaw_max)232 REAL, SAVE :: dens_rate233 !$OMP THREADPRIVATE(dens_rate)234 217 REAL :: wdens0 235 218 ! IM 080208 … … 239 222 REAL, DIMENSION (klon, klev) :: deltatw0 240 223 REAL, DIMENSION (klon, klev) :: deltaqw0 241 REAL, DIMENSION (klon, klev) :: te, qe 242 !! REAL, DIMENSION (klon) :: sigmaw1 224 REAL, DIMENSION (klon, klev) :: tb, qb 243 225 #ifdef ISO 244 226 REAL, DIMENSION (ntraciso,klon, klev) :: deltaxtw0 245 REAL, DIMENSION (ntraciso,klon, klev) :: xt e246 #endif 247 248 ! Variables liees a la dynamique de population 227 REAL, DIMENSION (ntraciso,klon, klev) :: xtb 228 #endif 229 230 ! Variables liees a la dynamique de population 1 249 231 REAL, DIMENSION(klon) :: act 250 232 REAL, DIMENSION(klon) :: rad_wk, tau_wk_inv 251 233 REAL, DIMENSION(klon) :: f_shear 252 234 REAL, DIMENSION(klon) :: drdt 253 REAL, DIMENSION(klon) :: d_sig_gen, d_sig_death, d_sig_col 235 236 ! Variables liees a la dynamique de population 2 237 REAL, DIMENSION(klon) :: cont_fact 238 239 ! Variables liees a la dynamique de population 3 240 REAL, DIMENSION(klon) :: arad_wk, irad_wk 241 242 !! REAL, DIMENSION(klon) :: d_sig_gen, d_sig_death, d_sig_col 254 243 REAL, DIMENSION(klon) :: wape1_act, wape2_act 255 244 LOGICAL, DIMENSION (klon) :: kill_wake 256 INTEGER, SAVE :: iflag_wk_act257 !$OMP THREADPRIVATE(iflag_wk_act)258 245 REAL :: drdt_pos 259 246 REAL :: tau_wk_inv_min 247 ! Some components of the tendencies of state variables 248 REAL, DIMENSION (klon) :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2 249 REAL, DIMENSION (klon) :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2 250 REAL, DIMENSION (klon) :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2 251 REAL, DIMENSION (klon) :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2 260 252 261 253 ! Variables pour les GW … … 266 258 267 259 ! Variables liees au calcul de hw 268 REAL, DIMENSION (klon) :: ptop _provis, ptop, ptop_new260 REAL, DIMENSION (klon) :: ptop 269 261 REAL, DIMENSION (klon) :: sum_dth 270 262 REAL, DIMENSION (klon) :: dthmin … … 278 270 ! Sub-timestep tendencies and related variables 279 271 REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw 280 REAL, DIMENSION (klon, klev) :: d_te, d_qe 281 REAL, DIMENSION (klon) :: d_awdens, d_wdens 282 REAL, DIMENSION (klon) :: d_sigmaw, alpha 272 REAL, DIMENSION (klon, klev) :: d_tb, d_qb 273 REAL, DIMENSION (klon) :: d_wdens, d_awdens, d_sigmaw, d_asigmaw 274 REAL, DIMENSION (klon) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd 275 REAL, DIMENSION (klon) :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd 276 REAL, DIMENSION (klon) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd 277 REAL, DIMENSION (klon) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd 278 REAL, DIMENSION (klon) :: agfl !! gust front length of active wakes 279 !! per unit area 280 REAL, DIMENSION (klon) :: alpha, alpha_tot 283 281 REAL, DIMENSION (klon) :: q0_min, q1_min 284 282 LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw 285 283 #ifdef ISO 286 284 REAL, DIMENSION (ntraciso,klon, klev) :: d_deltaxtw 287 REAL, DIMENSION (ntraciso,klon, klev) :: d_xte 288 #endif 289 REAL, SAVE :: epsilon 290 !$OMP THREADPRIVATE(epsilon) 291 DATA epsilon/1.E-15/ 285 REAL, DIMENSION (ntraciso,klon, klev) :: d_xtb 286 #endif 292 287 293 288 ! Autres variables internes 294 INTEGER ::isubstep, k, i 295 289 INTEGER ::isubstep, k, i, igout 290 291 REAL :: wdensmin 292 293 REAL :: sigmaw_targ 296 294 REAL :: wdens_targ 297 REAL :: sigmaw_targ 298 299 REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu 300 REAL, DIMENSION (klon) :: sum_dq, sum_rho 295 REAL :: d_sigmaw_targ 296 REAL :: d_wdens_targ 297 298 REAL, DIMENSION (klon) :: sum_thx, sum_tx, sum_qx, sum_thvx 299 REAL, DIMENSION (klon) :: sum_dq 301 300 REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn 302 REAL, DIMENSION (klon) :: av_th u, av_tu, av_qu, av_thvu303 REAL, DIMENSION (klon) :: av_dth, av_dq , av_rho301 REAL, DIMENSION (klon) :: av_thx, av_tx, av_qx, av_thvx 302 REAL, DIMENSION (klon) :: av_dth, av_dq 304 303 REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn 305 ! pas besoin d'isos ici 306 307 REAL, DIMENSION (klon, klev) :: rho, rhow 304 305 REAL, DIMENSION (klon, klev) :: rho 308 306 REAL, DIMENSION (klon, klev+1) :: rhoh 309 REAL, DIMENSION (klon, klev) :: rhow_moyen310 307 REAL, DIMENSION (klon, klev) :: zh 311 308 REAL, DIMENSION (klon, klev+1) :: zhh 312 REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2 313 314 REAL, DIMENSION (klon, klev) :: the, thu 309 310 REAL, DIMENSION (klon, klev) :: thb, thx 315 311 316 312 REAL, DIMENSION (klon, klev) :: omgbw … … 347 343 REAL, DIMENSION (klon, klev) :: detr 348 344 349 REAL, DIMENSION(klon) :: sigmaw_in ! pour les prints 350 REAL, DIMENSION(klon) :: awdens_in, wdens_in ! pour les prints 345 REAL, DIMENSION(klon) :: sigmaw_in, asigmaw_in ! pour les prints 346 REAL, DIMENSION(klon) :: wdens_in, awdens_in ! pour les prints 347 348 !!! LOGICAL :: phys_sub=.true. 349 LOGICAL :: phys_sub=.false. 350 351 LOGICAL :: first_call=.true. 352 353 354 !!-- variables liees au nouveau calcul de ptop et hw 355 REAL, DIMENSION (klon, klev) :: int_dth 356 REAL, DIMENSION (klon, klev) :: zzz, dzzz 357 REAL :: epsil 358 REAL, DIMENSION (klon) :: ptop1 359 INTEGER, DIMENSION (klon) :: ktop1 360 REAL, DIMENSION (klon) :: omega 361 REAL, DIMENSION (klon) :: h_zzz 362 363 !print*,'WAKE LJYFz' 351 364 352 365 ! ------------------------------------------------------------------------- 353 366 ! Initialisations 354 367 ! ------------------------------------------------------------------------- 355 356 ! print*, 'wake initialisations'357 !#ifdef ISOVERIF358 ! write(*,*) 'wake 358: entree'359 !#endif360 361 ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.362 ! -------------------------------------------------------------------------363 364 !! DATA wapecut, sigmad, hwmin/5., .02, 10./365 !! DATA wapecut, sigmad, hwmin/1., .02, 10./366 DATA sigmad, hwmin/.02, 10./367 !! DATA wdensmin/1.e-12/368 DATA wdensmin/1.e-14/369 ! cc nrlmd370 DATA sigmaw_max/0.4/371 DATA dens_rate/0.1/372 ! cc373 ! Longueur de maille (en m)374 ! -------------------------------------------------------------------------375 376 368 ! ALON = 3.e5 377 369 ! alon = 1.E6 … … 380 372 f_shear(:) = 1. ! 0. for strong shear, 1. for weak shear 381 373 382 383 374 ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) 384 375 385 ! coefgw : Coefficient pour les ondes de gravit é376 ! coefgw : Coefficient pour les ondes de gravite 386 377 ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) 387 ! wdens : Densit ésurfacique de poche froide378 ! wdens : Densite surfacique de poche froide 388 379 ! ------------------------------------------------------------------------- 389 380 … … 398 389 ! alpk = 0.5 399 390 ! alpk = 0.05 400 401 if (first) then402 403 igout = klon/2+1/klon404 405 crep_upper = 0.9406 crep_sol = 1.0407 408 ! Get wapecut from parameter file409 wapecut = 1.410 CALL getin_p('wapecut', wapecut)411 412 ! cc nrlmd Lecture du fichier wake_param.data413 stark=0.33414 CALL getin_p('stark',stark)415 cstart = stark*sqrt(2.*wapecut)416 417 alpk=0.25418 CALL getin_p('alpk',alpk)419 !jyg<420 !! wdens_ref=8.E-12421 !! CALL getin_p('wdens_ref',wdens_ref)422 wdens_ref(1)=8.E-12423 wdens_ref(2)=8.E-12424 CALL getin_p('wdens_ref_o',wdens_ref(1)) !wake number per unit area ; ocean425 CALL getin_p('wdens_ref_l',wdens_ref(2)) !wake number per unit area ; land426 !>jyg427 391 ! 428 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 429 !!!!!!!!! Population dynamics parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!! 430 !------------------------------------------------------------------------ 431 432 iflag_wk_pop_dyn = 0 433 CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed 434 ! and wdens prognostic 435 iflag_wk_act = 0 436 CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0. 437 ! 1: act(:)=1. 438 ! 2: act(:)=f(Wape) 439 440 rzero = 5000. 441 CALL getin_p('rzero_wk', rzero) 442 aa0 = 3.14*rzero*rzero 392 igout = klon/2+1/klon 443 393 ! 444 tau_cv = 4000. 445 CALL getin_p('tau_cv', tau_cv) 446 447 !------------------------------------------------------------------------ 448 449 coefgw=4. 450 CALL getin_p('coefgw',coefgw) 451 452 WRITE(*,*) 'stark=', stark 453 WRITE(*,*) 'alpk=', alpk 454 !jyg< 455 !! WRITE(*,*) 'wdens_ref=', wdens_ref 456 WRITE(*,*) 'wdens_ref_o=', wdens_ref(1) 457 WRITE(*,*) 'wdens_ref_l=', wdens_ref(2) 458 !>jyg 459 WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn 460 WRITE(*,*) 'iflag_wk_act',iflag_wk_act 461 WRITE(*,*) 'coefgw=', coefgw 462 463 flag_wk_check_trgl=.false. 464 CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl) 465 WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl 466 WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot' 467 iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1 468 CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl) 469 WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl 470 471 first=.false. 472 endif 394 ! sub-time-stepping parameters 395 dtimesub = dtime/wk_nsub 396 ! 397 IF (first_call) THEN 398 !!#define IOPHYS_WK 399 #undef IOPHYS_WK 400 #ifdef IOPHYS_WK 401 IF (phys_sub) THEN 402 call iophys_ini(dtimesub) 403 ELSE 404 call iophys_ini(dtime) 405 ENDIF 406 #endif 407 first_call = .false. 408 ENDIF !(first_call) 473 409 474 410 IF (iflag_wk_pop_dyn == 0) THEN … … 483 419 !>jyg 484 420 ENDIF ! (iflag_wk_pop_dyn == 0) 421 ! 422 IF (iflag_wk_pop_dyn >=1) THEN 423 IF (iflag_wk_pop_dyn == 3) THEN 424 wdensmin = wdensthreshold 425 ELSE 426 wdensmin = wdensinit 427 ENDIF 428 ENDIF 485 429 486 430 ! print*,'stark',stark … … 492 436 ! ------------------------------------------------------------------------- 493 437 494 delta_t_min = 0.2438 delta_t_min = 0.2 495 439 496 440 ! 1. - Save initial values, initialize tendencies, initialize output fields … … 503 447 !! deltatw0(i, k) = deltatw(i, k) 504 448 !! deltaqw0(i, k) = deltaqw(i, k) 505 !! t e(i, k) = te0(i, k)506 !! q e(i, k) = qe0(i, k)449 !! tb(i, k) = tb0(i, k) 450 !! qb(i, k) = qb0(i, k) 507 451 !! dtls(i, k) = 0. 508 452 !! dqls(i, k) = 0. 509 453 !! d_deltat_gw(i, k) = 0. 510 !! d_t e(i, k) = 0.511 !! d_q e(i, k) = 0.454 !! d_tb(i, k) = 0. 455 !! d_qb(i, k) = 0. 512 456 !! d_deltatw(i, k) = 0. 513 457 !! d_deltaqw(i, k) = 0. … … 521 465 deltatw0(:,:) = deltatw(:,:) 522 466 deltaqw0(:,:) = deltaqw(:,:) 523 t e(:,:) = te0(:,:)524 q e(:,:) = qe0(:,:)467 tb(:,:) = tb0(:,:) 468 qb(:,:) = qb0(:,:) 525 469 dtls(:,:) = 0. 526 470 dqls(:,:) = 0. 527 471 d_deltat_gw(:,:) = 0. 528 d_t e(:,:) = 0.529 d_q e(:,:) = 0.472 d_tb(:,:) = 0. 473 d_qb(:,:) = 0. 530 474 d_deltatw(:,:) = 0. 531 475 d_deltaqw(:,:) = 0. 532 476 d_deltatw2(:,:) = 0. 533 d_deltaqw2(:,:) = 0. 477 d_deltaqw2(:,:) = 0. 534 478 #ifdef ISO 535 479 deltaxtw0(:,:,:)= deltaxtw(:,:,:) 536 xt e(:,:,:) = xte0(:,:,:)480 xtb(:,:,:) = xtb0(:,:,:) 537 481 dxtls(:,:,:) = 0. 538 d_xt e(:,:,:) = 0.482 d_xtb(:,:,:) = 0. 539 483 d_deltaxtw(:,:,:) = 0. 540 484 d_deltaxtw2(:,:,:)=0. 541 xt1(:,:,:) = 0. 542 xt2(:,:,:)=0. 543 ! init non indispensable mais facilite verif iso 544 q1(:,:)=0.0 545 q2(:,:)=0.0 546 #endif 485 #endif 486 487 d_sig_gen2(:) = 0. 488 d_sig_death2(:) = 0. 489 d_sig_col2(:) = 0. 490 d_sig_spread2(:)= 0. 491 d_asig_death2(:) = 0. 492 d_asig_iicol2(:) = 0. 493 d_asig_aicol2(:) = 0. 494 d_asig_spread2(:)= 0. 495 d_asig_bnd2(:) = 0. 496 d_asigmaw2(:) = 0. 497 ! 498 d_dens_gen2(:) = 0. 499 d_dens_death2(:) = 0. 500 d_dens_col2(:) = 0. 501 d_dens_bnd2(:) = 0. 502 d_wdens2(:) = 0. 503 d_adens_bnd2(:) = 0. 504 d_awdens2(:) = 0. 505 d_adens_death2(:) = 0. 506 d_adens_icol2(:) = 0. 507 d_adens_acol2(:) = 0. 547 508 548 509 IF (iflag_wk_act == 0) THEN … … 555 516 !! sigmaw_in(i) = sigmaw(i) 556 517 !! END DO 557 sigmaw_in(:) = sigmaw(:) 518 sigmaw_in(:) = sigmaw(:) 519 asigmaw_in(:) = asigmaw(:) 558 520 !>jyg 559 560 ! sigmaw1=sigmaw561 ! IF (sigd_con.GT.sigmaw1) THEN562 ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con563 ! ENDIF564 IF (iflag_wk_pop_dyn >=1) THEN565 DO i = 1, klon566 wdens_targ = max(wdens(i),wdensmin)567 d_wdens2(i) = wdens_targ - wdens(i)568 wdens(i) = wdens_targ569 END DO570 ELSE571 DO i = 1, klon572 d_awdens2(i) = 0.573 d_wdens2(i) = 0.574 END DO575 ENDIF ! (iflag_wk_pop_dyn >=1)576 !577 DO i = 1, klon578 ! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i))579 !jyg<580 !! sigmaw(i) = amax1(sigmaw(i), sigmad)581 !! sigmaw(i) = amin1(sigmaw(i), 0.99)582 sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)583 d_sigmaw2(i) = sigmaw_targ - sigmaw(i)584 sigmaw(i) = sigmaw_targ585 !>jyg586 END DO587 588 521 ! 589 522 IF (iflag_wk_pop_dyn >= 1) THEN … … 595 528 ENDIF ! (iflag_wk_pop_dyn >= 1) 596 529 530 531 ! sigmaw1=sigmaw 532 ! IF (sigd_con.GT.sigmaw1) THEN 533 ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con 534 ! ENDIF 535 IF (iflag_wk_pop_dyn >= 1) THEN 536 DO i = 1, klon 537 d_dens_gen2(i) = 0. 538 d_dens_death2(i) = 0. 539 d_dens_col2(i) = 0. 540 d_awdens2(i) = 0. 541 IF (wdens(i) < wdensthreshold) THEN 542 !! wdens_targ = max(wdens(i),wdensmin) 543 wdens_targ = max(wdens(i),wdensinit) 544 d_dens_bnd2(i) = wdens_targ - wdens(i) 545 d_wdens2(i) = wdens_targ - wdens(i) 546 wdens(i) = wdens_targ 547 ELSE 548 d_dens_bnd2(i) = 0. 549 d_wdens2(i) = 0. 550 ENDIF !! (wdens(i) < wdensthreshold) 551 END DO 552 IF (iflag_wk_pop_dyn >= 2) THEN 553 DO i = 1, klon 554 IF (awdens(i) < wdensthreshold) THEN 555 !! wdens_targ = min(max(awdens(i),wdensmin),wdens(i)) 556 wdens_targ = min(max(awdens(i),wdensinit),wdens(i)) 557 d_adens_bnd2(i) = wdens_targ - awdens(i) 558 d_awdens2(i) = wdens_targ - awdens(i) 559 awdens(i) = wdens_targ 560 ELSE 561 wdens_targ = min(awdens(i), wdens(i)) 562 d_adens_bnd2(i) = wdens_targ - awdens(i) 563 d_awdens2(i) = wdens_targ - awdens(i) 564 awdens(i) = wdens_targ 565 ENDIF 566 END DO 567 ENDIF ! (iflag_wk_pop_dyn >= 2) 568 ELSE 569 DO i = 1, klon 570 d_awdens2(i) = 0. 571 d_wdens2(i) = 0. 572 END DO 573 ENDIF ! (iflag_wk_pop_dyn >= 1) 574 ! 575 DO i = 1, klon 576 sigmaw_targ = min(max(sigmaw(i), sigmad),0.99) 577 d_sig_bnd2(i) = sigmaw_targ - sigmaw(i) 578 d_sigmaw2(i) = sigmaw_targ - sigmaw(i) 579 sigmaw(i) = sigmaw_targ 580 END DO 581 ! 582 IF (iflag_wk_pop_dyn == 3) THEN 583 DO i = 1, klon 584 IF ((wdens(i)-awdens(i)) <= smallestreal) THEN 585 sigmaw_targ = sigmaw(i) 586 ELSE 587 sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i)) 588 ENDIF 589 d_asig_bnd2(i) = sigmaw_targ - asigmaw(i) 590 d_asigmaw2(i) = sigmaw_targ - asigmaw(i) 591 asigmaw(i) = sigmaw_targ 592 END DO 593 ENDIF ! (iflag_wk_pop_dyn == 3) 594 597 595 wape(:) = 0. 598 596 wape2(:) = 0. 599 597 d_sigmaw(:) = 0. 598 d_asigmaw(:) = 0. 600 599 ktopw(:) = 0 601 600 ! 602 601 !<jyg 603 602 dth(:,:) = 0. 604 t u(:,:) = 0.605 q u(:,:) = 0.603 tx(:,:) = 0. 604 qx(:,:) = 0. 606 605 dtke(:,:) = 0. 607 606 dqke(:,:) = 0. 608 spread(:,:) = 0.607 wkspread(:,:) = 0. 609 608 omgbdth(:,:) = 0. 610 609 omg(:,:) = 0. … … 624 623 omgbdq(:,:) = 0. 625 624 #ifdef ISO 626 xt u(:,:,:) = 0.625 xtx(:,:,:) = 0. 627 626 dxtke(:,:,:) = 0. 628 627 omgbdxt(:,:,:) = 0. … … 653 652 ktop(i) = 0 654 653 kupper(i) = 0 655 sum_th u(i) = 0.656 sum_t u(i) = 0.657 sum_q u(i) = 0.658 sum_thv u(i) = 0.654 sum_thx(i) = 0. 655 sum_tx(i) = 0. 656 sum_qx(i) = 0. 657 sum_thvx(i) = 0. 659 658 sum_dth(i) = 0. 660 659 sum_dq(i) = 0. 661 sum_rho(i) = 0.662 660 sum_dtdwn(i) = 0. 663 661 sum_dqdwn(i) = 0. 664 662 665 av_th u(i) = 0.666 av_t u(i) = 0.667 av_q u(i) = 0.668 av_thv u(i) = 0.663 av_thx(i) = 0. 664 av_tx(i) = 0. 665 av_qx(i) = 0. 666 av_thvx(i) = 0. 669 667 av_dth(i) = 0. 670 668 av_dq(i) = 0. 671 av_rho(i) = 0.672 669 av_dtdwn(i) = 0. 673 670 av_dqdwn(i) = 0. 674 671 ! pas besoin d'isos ici 675 672 END DO 676 677 678 #ifdef ISOVERIF679 ! TODO680 #endif681 673 682 674 ! Distance between wakes … … 688 680 DO k = 1, klev 689 681 DO i = 1, klon 690 ! write(*,*)'wake 1',i,k, rd,te(i,k)691 rho(i, k) = p(i, k)/( rd*te(i,k))682 ! write(*,*)'wake 1',i,k,RD,tb(i,k) 683 rho(i, k) = p(i, k)/(RD*tb(i,k)) 692 684 ! write(*,*)'wake 2',rho(i,k) 693 685 IF (k==1) THEN 694 ! write(*,*)'wake 3',i,k,rd,t e(i,k)695 rhoh(i, k) = ph(i, k)/( rd*te(i,k))696 ! write(*,*)'wake 4',i,k,rd,t e(i,k)686 ! write(*,*)'wake 3',i,k,rd,tb(i,k) 687 rhoh(i, k) = ph(i, k)/(RD*tb(i,k)) 688 ! write(*,*)'wake 4',i,k,rd,tb(i,k) 697 689 zhh(i, k) = 0 698 690 ELSE 699 ! write(*,*)'wake 5',rd,(t e(i,k)+te(i,k-1))700 rhoh(i, k) = ph(i, k)*2./( rd*(te(i,k)+te(i,k-1)))691 ! write(*,*)'wake 5',rd,(tb(i,k)+tb(i,k-1)) 692 rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1))) 701 693 ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1) 702 zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)* rg) + zhh(i, k-1)694 zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1) 703 695 END IF 704 696 ! write(*,*)'wake 7',ppi(i,k) 705 the(i, k) = te(i, k)/ppi(i, k) 706 thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 707 tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) 708 qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) 709 ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k))) 710 rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) 697 thb(i, k) = tb(i, k)/ppi(i, k) 698 thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 699 tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i) 700 qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i) 701 ! write(*,*)'wake 8',(RD*(tb(i,k)+deltatw(i,k))) 711 702 dth(i, k) = deltatw(i, k)/ppi(i, k) 712 703 #ifdef ISO 713 704 do ixt=1,ntraciso 714 xt u(ixt,i,k) = xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)705 xtx(ixt,i,k) = xtb(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i) 715 706 enddo !do ixt=1,ntraciso 716 707 #ifdef ISOVERIF 717 708 if (iso_eau.gt.0) then 718 709 call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 723a') 719 call iso_verif_egalite(q u(i,k),xtu(iso_eau,i,k),'wake 723b')710 call iso_verif_egalite(qx(i,k),xtx(iso_eau,i,k),'wake 723b') 720 711 endif 721 712 if (iso_HDO.gt.0) then 722 if (iso_verif_aberrant_enc_choix_nostop(xt u(iso_hdo,i,k),qu(i,k),ridicule,deltalim, &723 'wake 723c xt u').eq.1) then713 if (iso_verif_aberrant_enc_choix_nostop(xtx(iso_hdo,i,k),qx(i,k),ridicule,deltalim, & 714 'wake 723c xtx').eq.1) then 724 715 stop 725 716 endif … … 730 721 END DO 731 722 732 733 734 723 DO k = 1, klev - 1 735 724 DO i = 1, klon … … 737 726 n2(i, k) = 0 738 727 ELSE 739 n2(i, k) = amax1(0., - rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ &728 n2(i, k) = amax1(0., -RG**2/thb(i,k)*rho(i,k)*(thb(i,k+1)-thb(i,k-1))/ & 740 729 (p(i,k+1)-p(i,k-1))) 741 730 END IF … … 754 743 END DO 755 744 756 ! Calcul de la masse volumique moyenne de la colonne (bdlmd) 757 ! ----------------------------------------------------------------- 758 759 DO k = 1, klev 760 DO i = 1, klon 761 epaisseur1(i, k) = 0. 762 epaisseur2(i, k) = 0. 763 END DO 764 END DO 765 766 DO i = 1, klon 767 epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. 768 epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. 769 rhow_moyen(i, 1) = rhow(i, 1) 770 END DO 771 772 DO k = 2, klev 773 DO i = 1, klon 774 epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1. 775 epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k) 776 rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* & 777 epaisseur1(i,k))/epaisseur2(i, k) 778 END DO 779 END DO 780 781 745 782 746 ! Choose an integration bound well above wake top 783 747 ! ----------------------------------------------------------------- 784 748 785 ! Pupper = 50000. ! melting level786 ! Pupper = 60000.787 ! Pupper = 80000. ! essais pour case_e788 DO i = 1, klon789 pupper(i) = 0.6*ph(i, 1)790 pupper(i) = max(pupper(i), 45000.)791 ! cc Pupper(i) = 60000.792 END DO793 794 795 749 ! Determine Wake top pressure (Ptop) from buoyancy integral 796 750 ! -------------------------------------------------------- 797 751 798 ! -1/ Pressure of the level where dth becomes less than delta_t_min. 799 800 DO i = 1, klon 801 ptop_provis(i) = ph(i, 1) 802 END DO 803 DO k = 2, klev 804 DO i = 1, klon 805 806 ! IM v3JYG; ptop_provis(i).LT. ph(i,1) 807 808 IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. & 809 ptop_provis(i)==ph(i,1)) THEN 810 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- & 811 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 812 END IF 813 END DO 814 END DO 815 816 ! -2/ dth integral 817 818 DO i = 1, klon 819 sum_dth(i) = 0. 820 dthmin(i) = -delta_t_min 821 z(i) = 0. 822 END DO 823 824 DO k = 1, klev 825 DO i = 1, klon 826 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) 827 IF (dz(i)>0) THEN 828 z(i) = z(i) + dz(i) 829 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 830 dthmin(i) = amin1(dthmin(i), dth(i,k)) 831 END IF 832 END DO 833 END DO 834 835 ! -3/ height of triangle with area= sum_dth and base = dthmin 836 837 DO i = 1, klon 838 hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) 839 hw0(i) = amax1(hwmin, hw0(i)) 840 END DO 841 842 ! -4/ now, get Ptop 843 844 DO i = 1, klon 845 z(i) = 0. 846 ptop(i) = ph(i, 1) 847 END DO 848 849 DO k = 1, klev 850 DO i = 1, klon 851 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i)) 852 IF (dz(i)>0) THEN 853 z(i) = z(i) + dz(i) 854 ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) 855 END IF 856 END DO 857 END DO 858 859 IF (prt_level>=10) THEN 860 PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout) 861 ENDIF 862 863 864 ! -5/ Determination de ktop et kupper 865 866 DO k = klev, 1, -1 867 DO i = 1, klon 868 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 869 IF (ph(i,k+1)<pupper(i)) kupper(i) = k 870 END DO 871 END DO 872 873 ! On evite kupper = 1 et kupper = klev 874 DO i = 1, klon 875 kupper(i) = max(kupper(i), 2) 876 kupper(i) = min(kupper(i), klev-1) 877 END DO 878 879 880 ! -6/ Correct ktop and ptop 881 882 DO i = 1, klon 883 ptop_new(i) = ptop(i) 884 END DO 885 DO k = klev, 2, -1 886 DO i = 1, klon 887 IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & 888 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 889 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & 890 k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 891 END IF 892 END DO 893 END DO 894 895 DO i = 1, klon 896 ptop(i) = ptop_new(i) 897 END DO 898 899 DO k = klev, 1, -1 900 DO i = 1, klon 901 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 902 END DO 903 END DO 904 905 IF (prt_level>=10) THEN 906 PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) 752 Do i=1, klon 753 wk_adv(i) = .True. 754 Enddo 755 Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, & 756 dth, hw0, rho, delta_t_min, & 757 ktop, wk_adv, h_zzz, ptop1, ktop1) 758 759 !!print'("pkupper APPEL ",7i6)',0,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper 760 761 IF (prt_level>=10) THEN 762 PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) 907 763 ENDIF 908 764 … … 915 771 deltaqw(i, k) = 0. 916 772 d_deltatw2(i,k) = -deltatw0(i,k) 773 d_deltaqw2(i,k) = -deltaqw0(i,k) 917 774 #ifdef ISO 918 775 do ixt=1,ntraciso … … 939 796 ! -------------------------------------- 940 797 941 ! Initialize sum_thv uto 1st level virt. pot. temp.798 ! Initialize sum_thvx to 1st level virt. pot. temp. 942 799 943 800 DO i = 1, klon 944 801 z(i) = 1. 945 802 dz(i) = 1. 946 sum_thv u(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)803 sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i) 947 804 sum_dth(i) = 0. 948 805 END DO … … 950 807 DO k = 1, klev 951 808 DO i = 1, klon 952 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)* rg)809 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG) 953 810 IF (dz(i)>0) THEN 811 ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau 954 812 z(i) = z(i) + dz(i) 955 sum_th u(i) = sum_thu(i) + thu(i, k)*dz(i)956 sum_t u(i) = sum_tu(i) + tu(i, k)*dz(i)957 sum_q u(i) = sum_qu(i) + qu(i, k)*dz(i)958 sum_thv u(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)813 sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i) 814 sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i) 815 sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i) 816 sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i) 959 817 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 960 818 sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) 961 sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)962 819 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) 963 820 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) … … 979 836 980 837 DO i = 1, klon 981 av_th u(i) = sum_thu(i)/hw0(i)982 av_t u(i) = sum_tu(i)/hw0(i)983 av_q u(i) = sum_qu(i)/hw0(i)984 av_thv u(i) = sum_thvu(i)/hw0(i)838 av_thx(i) = sum_thx(i)/hw0(i) 839 av_tx(i) = sum_tx(i)/hw0(i) 840 av_qx(i) = sum_qx(i)/hw0(i) 841 av_thvx(i) = sum_thvx(i)/hw0(i) 985 842 ! av_thve = sum_thve/hw0 986 843 av_dth(i) = sum_dth(i)/hw0(i) 987 844 av_dq(i) = sum_dq(i)/hw0(i) 988 av_rho(i) = sum_rho(i)/hw0(i)989 845 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 990 846 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 991 847 992 wape(i) = - rg*hw0(i)*(av_dth(i)+ &993 epsim1*(av_th u(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)848 wape(i) = -RG*hw0(i)*(av_dth(i)+ & 849 epsim1*(av_thx(i)*av_dq(i)+av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i) 994 850 995 851 END DO 852 #ifdef IOPHYS_WK 853 IF (.not.phys_sub) CALL iophys_ecrit('wape_a',1,'wape_a','J/kg',wape) 854 #endif 996 855 997 856 ! 2.2 Prognostic variable update … … 1020 879 DO i = 1, klon 1021 880 IF (wape(i)<0.) THEN 881 !! sigmaw(i) = amax1(sigmad, sigd_con(i)) 882 sigmaw_targ = max(sigmad, sigd_con(i)) 883 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) 884 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 885 sigmaw(i) = sigmaw_targ 886 ENDIF !! (wape(i)<0.) 887 ENDDO 888 ! 889 IF (iflag_wk_pop_dyn == 3) THEN 890 DO i = 1, klon 891 IF (wape(i)<0.) THEN 892 sigmaw_targ = max(sigmad, sigd_con(i)) 893 d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) 894 d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) 895 asigmaw(i) = sigmaw_targ 896 ENDIF !! (wape(i)<0.) 897 ENDDO 898 ENDIF !! (iflag_wk_pop_dyn == 3) 899 900 DO i = 1, klon 901 IF (wape(i)<0.) THEN 1022 902 wape(i) = 0. 1023 903 cstar(i) = 0. 1024 904 hw(i) = hwmin 1025 !jyg<1026 !! sigmaw(i) = amax1(sigmad, sigd_con(i))1027 sigmaw_targ = max(sigmad, sigd_con(i))1028 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)1029 sigmaw(i) = sigmaw_targ1030 !>jyg1031 905 fip(i) = 0. 1032 906 gwake(i) = .FALSE. … … 1037 911 END IF 1038 912 END DO 1039 913 ! 1040 914 1041 915 ! Check qx and qw positivity 1042 916 ! -------------------------- 1043 917 DO i = 1, klon 1044 q0_min(i) = min((q e(i,1)-sigmaw(i)*deltaqw(i,1)), &1045 (q e(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))918 q0_min(i) = min((qb(i,1)-sigmaw(i)*deltaqw(i,1)), & 919 (qb(i,1)+(1.-sigmaw(i))*deltaqw(i,1))) 1046 920 END DO 1047 921 DO k = 2, klev 1048 922 DO i = 1, klon 1049 q1_min(i) = min((q e(i,k)-sigmaw(i)*deltaqw(i,k)), &1050 (q e(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))923 q1_min(i) = min((qb(i,k)-sigmaw(i)*deltaqw(i,k)), & 924 (qb(i,k)+(1.-sigmaw(i))*deltaqw(i,k))) 1051 925 IF (q1_min(i)<=q0_min(i)) THEN 1052 926 q0_min(i) = q1_min(i) … … 1058 932 ok_qx_qw(i) = q0_min(i) >= 0. 1059 933 alpha(i) = 1. 934 alpha_tot(i) = 1. 1060 935 END DO 1061 936 … … 1070 945 ! ----------------- 1071 946 1072 nsub = 10 1073 dtimesub = dtime/nsub 1074 1075 ! ------------------------------------------------------------ 1076 DO isubstep = 1, nsub 1077 ! ------------------------------------------------------------ 1078 947 ! wk_nsub and dtimesub definitions moved to begining of routine. 948 !! wk_nsub = 10 949 !! dtimesub = dtime/wk_nsub 950 951 952 ! ------------------------------------------------------------------------ 953 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 954 ! ------------------------------------------------------------------------ 955 ! 956 DO isubstep = 1, wk_nsub 957 ! 958 ! ------------------------------------------------------------------------ 1079 959 ! wk_adv is the logical flag enabling wake evolution in the time advance 1080 960 ! loop … … 1085 965 PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', & 1086 966 isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) 967 1087 968 ENDIF 1088 969 1089 970 ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement 1090 ! n égatif de ktop àkupper --------1091 ! cc On calcule pour cela une densit éwdens0 pour laquelle on971 ! negatif de ktop a kupper -------- 972 ! cc On calcule pour cela une densite wdens0 pour laquelle on 1092 973 ! aurait un entrainement nul --- 1093 974 !jyg< … … 1096 977 ! des descentes unsaturees. Nous faisons alors l'hypothese que la 1097 978 ! convection profonde cree directement de nouvelles poches, sans passer 1098 ! par les thermiques. La nouvelle valeur de wdens est alors impos ée.979 ! par les thermiques. La nouvelle valeur de wdens est alors imposee. 1099 980 1100 981 DO i = 1, klon … … 1102 983 ! c $ isubstep,wk_adv(i),cstar(i),wape(i) 1103 984 IF (wk_adv(i) .AND. cstar(i)>0.01) THEN 1104 omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & 1105 rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) 985 IF ( iflag_wk_profile == 0 ) THEN 986 omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + & 987 RG*amup(i, kupper(i)+1)/(1.-sigmaw(i)) 988 ELSE 989 omg(i, kupper(i)+1)=0. 990 ENDIF 1106 991 wdens0 = (sigmaw(i)/(4.*3.14))* & 1107 992 ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2) … … 1112 997 IF (wdens(i)<=wdens0*1.1) THEN 1113 998 IF (iflag_wk_pop_dyn >= 1) THEN 999 d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i) 1114 1000 d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i) 1115 1001 ENDIF … … 1119 1005 END DO 1120 1006 1007 IF (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) THEN 1008 !!-------------------------------------------------------- 1009 !!Bug : computing gfl and rad_wk before changing sigmaw 1010 !! This bug exists only for iflag_wk_pop_dyn=0. Otherwise, gfl and rad_wk 1011 !! are computed within wake_popdyn 1012 !!-------------------------------------------------------- 1013 DO i = 1, klon 1014 IF (wk_adv(i)) THEN 1015 gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) 1016 rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i))) 1017 END IF 1018 END DO 1019 ENDIF ! (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) 1020 !!-------------------------------------------------------- 1021 1121 1022 DO i = 1, klon 1122 1023 IF (wk_adv(i)) THEN 1123 gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))1124 rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))1125 !jyg<1126 !! sigmaw(i) = amin1(sigmaw(i), sigmaw_max)1127 1024 sigmaw_targ = min(sigmaw(i), sigmaw_max) 1128 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 1025 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) 1026 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 1129 1027 sigmaw(i) = sigmaw_targ 1130 !>jyg 1131 END IF 1132 END DO 1028 END IF 1029 END DO 1030 1031 IF (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) THEN 1032 !!-------------------------------------------------------- 1033 !!Fix : computing gfl and rad_wk after changing sigmaw 1034 !!-------------------------------------------------------- 1035 DO i = 1, klon 1036 IF (wk_adv(i)) THEN 1037 gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) 1038 rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i))) 1039 END IF 1040 END DO 1041 ENDIF ! (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) 1042 !!-------------------------------------------------------- 1133 1043 1134 1044 IF (iflag_wk_pop_dyn >= 1) THEN 1135 !The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.1136 !Here, it has to be set to zero.1045 ! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0. 1046 ! Here, it has to be set to zero. 1137 1047 death_rate(:) = 0. 1138 1139 IF (iflag_wk_act ==2) THEN 1048 ENDIF 1049 1050 IF (iflag_wk_pop_dyn >= 3) THEN 1140 1051 DO i = 1, klon 1141 1052 IF (wk_adv(i)) THEN 1142 wape1_act(i) = abs(cin(i)) 1143 wape2_act(i) = 2.*wape1_act(i) + 1. 1144 act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) )) 1145 ENDIF ! (wk_adv(i)) 1146 ENDDO 1147 ENDIF ! (iflag_wk_act ==2) 1148 1149 1150 DO i = 1, klon 1151 IF (wk_adv(i)) THEN 1152 !! tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.) 1153 tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.) 1154 tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub) 1155 drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / & 1156 (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i)) 1157 !! (1 - 2*sigmaw(i)*(1.-f_shear(i))) 1158 drdt_pos=max(drdt(i),0.) 1159 1160 !! d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) & 1161 !! - wdens(i)*tau_wk_inv_min & 1162 !! - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub 1163 d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub 1164 d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min - & 1165 2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub 1166 d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i)) 1167 1168 !! d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) & 1169 !! + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) & 1170 !! - sigmaw(i)*tau_wk_inv_min )*dtimesub 1171 d_sig_gen(i) = wgen(i)*aa0 1172 d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min 1173 !! d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos 1174 d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos 1175 d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub 1176 d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) 1053 sigmaw_targ = min(asigmaw(i), sigmaw_max) 1054 d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) 1055 d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) 1056 asigmaw(i) = sigmaw_targ 1177 1057 ENDIF 1178 1058 ENDDO 1179 1180 IF (prt_level >= 10) THEN 1181 print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', & 1182 cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) 1183 print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', & 1184 wdens(1), awdens(1), act(1), d_awdens(1) 1185 print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', & 1186 wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1) 1187 print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & 1188 d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) 1189 ENDIF 1059 ENDIF 1060 1061 !!-------------------------------------------------------- 1062 !!-------------------------------------------------------- 1063 IF (iflag_wk_pop_dyn == 1) THEN 1064 ! 1065 CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, & 1066 wdensmin, & 1067 dtimesub, gfl, rad_wk, f_shear, drdt_pos, & 1068 d_awdens, d_wdens, d_sigmaw, & 1069 iflag_wk_act, wk_adv, cin, wape, & 1070 drdt, & 1071 d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & 1072 d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & 1073 d_wdens_targ, d_sigmaw_targ) 1074 1190 1075 1191 ELSE ! (iflag_wk_pop_dyn >= 1) 1192 1076 !!-------------------------------------------------------- 1077 ELSEIF (iflag_wk_pop_dyn == 2) THEN 1078 ! 1079 CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, & 1080 wdensmin, & 1081 sigmaw, wdens, awdens, & !! state variables 1082 gfl, cstar, cin, wape, rad_wk, & 1083 d_sigmaw, d_wdens, d_awdens, & !! tendencies 1084 cont_fact, & 1085 d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & 1086 d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & 1087 d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd ) 1088 sigmaw=sigmaw-d_sigmaw 1089 wdens=wdens-d_wdens 1090 awdens=awdens-d_awdens 1091 1092 !!-------------------------------------------------------- 1093 ELSEIF (iflag_wk_pop_dyn == 3) THEN 1094 #ifdef IOPHYS_WK 1095 IF (phys_sub) THEN 1096 CALL iophys_ecrit('ptop',1,'ptop','Pa',ptop) 1097 CALL iophys_ecrit('sigmaw',1,'sigmaw','',sigmaw) 1098 CALL iophys_ecrit('asigmaw',1,'asigmaw','',asigmaw) 1099 CALL iophys_ecrit('wdens',1,'wdens','1/m2',wdens) 1100 CALL iophys_ecrit('awdens',1,'awdens','1/m2',awdens) 1101 CALL iophys_ecrit('rad_wk',1,'rad_wk','m',rad_wk) 1102 CALL iophys_ecrit('arad_wk',1,'arad_wk','m',arad_wk) 1103 CALL iophys_ecrit('irad_wk',1,'irad_wk','m',irad_wk) 1104 ENDIF 1105 #endif 1106 ! 1107 CALL wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, & 1108 wdensmin, & 1109 sigmaw, asigmaw, wdens, awdens, & !! state variables 1110 gfl, agfl, cstar, cin, wape, & 1111 rad_wk, arad_wk, irad_wk, & 1112 d_sigmaw, d_asigmaw, d_wdens, d_awdens, & !! tendencies 1113 d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & 1114 d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, & 1115 d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & 1116 d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd ) 1117 sigmaw=sigmaw-d_sigmaw 1118 asigmaw=asigmaw-d_asigmaw 1119 wdens=wdens-d_wdens 1120 awdens=awdens-d_awdens 1121 1122 !!-------------------------------------------------------- 1123 ELSEIF (iflag_wk_pop_dyn == 0) THEN 1124 1193 1125 ! cc nrlmd 1194 1126 1195 1127 DO i = 1, klon 1196 1128 IF (wk_adv(i)) THEN 1197 ! cc nrlmd Introduction du taux de mortalité des poches et 1129 1130 ! cc nrlmd Introduction du taux de mortalite des poches et 1198 1131 ! test sur sigmaw_max=0.4 1199 1132 ! cc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub … … 1220 1153 END DO 1221 1154 1222 ENDIF ! (iflag_wk_pop_dyn >= 1) 1223 1224 1155 ENDIF ! (iflag_wk_pop_dyn == 1) ... ELSEIF (iflag_wk_pop_dyn == 0) 1156 !!-------------------------------------------------------- 1157 !!-------------------------------------------------------- 1158 1159 #ifdef IOPHYS_WK 1160 IF (phys_sub) THEN 1161 CALL iophys_ecrit('wdensa',1,'wdensa','m',wdens) 1162 CALL iophys_ecrit('awdensa',1,'awdensa','m',awdens) 1163 CALL iophys_ecrit('sigmawa',1,'sigmawa','m',sigmaw) 1164 CALL iophys_ecrit('asigmawa',1,'asigmawa','m',asigmaw) 1165 ENDIF 1166 #endif 1225 1167 ! calcul de la difference de vitesse verticale poche - zone non perturbee 1226 1168 ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg 1227 1169 ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit 1228 ! IM 060208 au niveau k=1.. ?1170 ! IM 060208 au niveau k=1... 1229 1171 !JYG 161013 Correction : maintenant omg est dimensionne a klev. 1230 1172 DO k = 1, klev … … 1254 1196 DO i = 1, klon 1255 1197 IF (wk_adv(i) .AND. k<=ktop(i)) THEN 1256 dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)* rg)1198 dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*RG) 1257 1199 z(i) = z(i) + dz(i) 1258 1200 dp_deltomg(i, k) = dp_deltomg(i, 1) … … 1264 1206 DO i = 1, klon 1265 1207 IF (wk_adv(i)) THEN 1266 dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))* rg)1208 dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*RG) 1267 1209 ztop(i) = z(i) + dztop(i) 1268 1210 omgtop(i) = dp_deltomg(i, 1)*ztop(i) … … 1282 1224 DO i = 1, klon 1283 1225 IF (wk_adv(i)) THEN 1284 omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i) 1285 dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) 1226 omgtop(i) = -rho(i, ktop(i))*RG*omgtop(i) 1227 !! LJYF dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) 1228 dp_deltomg(i, 1) = omgtop(i)/min(ptop(i)-ph(i,1),-smallestreal) 1286 1229 END IF 1287 1230 END DO … … 1290 1233 DO i = 1, klon 1291 1234 IF (wk_adv(i) .AND. k<=ktop(i)) THEN 1292 omg(i, k) = -rho(i, k)* rg*omg(i, k)1235 omg(i, k) = -rho(i, k)*RG*omg(i, k) 1293 1236 dp_deltomg(i, k) = dp_deltomg(i, 1) 1294 1237 END IF … … 1300 1243 DO i = 1, klon 1301 1244 IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN 1302 omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & 1303 rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) 1245 IF ( iflag_wk_profile == 0 ) THEN 1246 omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + & 1247 RG*amup(i, kupper(i)+1)/(1.-sigmaw(i)) 1248 ELSE 1249 omg(i, kupper(i)+1) = 0. 1250 ENDIF 1304 1251 dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ & 1305 1252 (ptop(i)-pupper(i)) … … 1308 1255 1309 1256 ! c DO i=1,klon 1310 ! c print*,'Pente entre 0 et kupper (r éférence)'1257 ! c print*,'Pente entre 0 et kupper (reference)' 1311 1258 ! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) 1312 1259 ! c print*,'Pente entre ktop et kupper' … … 1384 1331 IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN 1385 1332 dth(i, k) = deltatw(i, k)/ppi(i, k) 1386 th1(i, k) = th e(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area1387 th2(i, k) = th e(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake1388 q1(i, k) = q e(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area1389 q2(i, k) = q e(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake1333 th1(i, k) = thb(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area 1334 th2(i, k) = thb(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake 1335 q1(i, k) = qb(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area 1336 q2(i, k) = qb(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake 1390 1337 #ifdef ISO 1391 1338 do ixt=1,ntraciso 1392 xt1(ixt,i,k) = xt e(ixt,i,k) -sigmaw(i) *deltaxtw(ixt,i,k) ! undisturbed area1393 xt2(ixt,i,k) = xt e(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake1339 xt1(ixt,i,k) = xtb(ixt,i,k) -sigmaw(i) *deltaxtw(ixt,i,k) ! undisturbed area 1340 xt2(ixt,i,k) = xtb(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake 1394 1341 enddo 1395 1342 #endif … … 1469 1416 END DO 1470 1417 1471 IF (prt_level>=10) THEN 1472 PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,klev) 1473 PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,klev) 1474 PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,klev) 1475 PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,klev) 1418 !! IF (prt_level>=10) THEN 1419 IF (prt_level>=10 .and. wk_adv(igout)) THEN 1420 PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout)) 1421 PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout)) 1422 PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout)) 1423 PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout)) 1476 1424 ENDIF 1477 1425 … … 1489 1437 (1.-alpha_up(i,k))*omgbdth(i,k)- & 1490 1438 alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k) 1491 ! print*,'d_d eltatw=', k, d_deltatw(i,k)1439 ! print*,'d_d,k_ptop_provis(i)eltatw=', k, d_deltatw(i,k) 1492 1440 1493 1441 d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & … … 1523 1471 ! C 1524 1472 ! ----------------------------------------------------------------- 1525 d_t e(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &1473 d_tb(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- & 1526 1474 rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/ & 1527 1475 (ph(i,k)-ph(i,k+1)) & … … 1529 1477 (ph(i,k)-ph(i,k+1)) )*ppi(i, k) 1530 1478 1531 d_q e(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- &1479 d_qb(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- & 1532 1480 rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ & 1533 1481 (ph(i,k)-ph(i,k+1)) & … … 1536 1484 #ifdef ISO 1537 1485 do ixt=1,ntraciso 1538 d_xt e(ixt,i,k) = dtimesub*( &1486 d_xtb(ixt,i,k) = dtimesub*( & 1539 1487 ( rre1(i)*omg(i,k )*sigmaw(i) *d_xt1(ixt,i,k) & 1540 1488 -rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1) ) & … … 1546 1494 #endif 1547 1495 ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN 1548 d_t e(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)1549 1550 d_q e(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))1496 d_tb(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k) 1497 1498 d_qb(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1))) 1551 1499 1552 1500 #ifdef ISO 1553 1501 do ixt=1,ntraciso 1554 d_xt e(ixt,i,k) = dtimesub*( &1502 d_xtb(ixt,i,k) = dtimesub*( & 1555 1503 ( rre1(i)*omg(i,k )*sigmaw(i) *d_xt1(ixt,i,k) & 1556 1504 /(Ph(i,k)-Ph(i,k+1))) & … … 1602 1550 1603 1551 1604 ! Coefficient de r épartition1552 ! Coefficient de repartition 1605 1553 1606 1554 crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ & 1607 1555 (ph(i,kupper(i))-ph(i,1)) 1608 1556 crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ & 1609 (p (i,1)-ph(i,kupper(i)))1557 (ph(i,1)-ph(i,kupper(i))) 1610 1558 1611 1559 … … 1646 1594 ! 1647 1595 1648 ! cc nrlmd Prise en compte du taux de mortalit é1649 ! cc D éfinitions de entr, detr1596 ! cc nrlmd Prise en compte du taux de mortalite 1597 ! cc Definitions de entr, detr 1650 1598 !jyg< 1651 1599 !! detr(i, k) = 0. … … 1657 1605 sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) 1658 1606 !>jyg 1659 spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)1660 1661 ! cc spread(i,k) =1607 wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i) 1608 1609 ! cc wkspread(i,k) = 1662 1610 ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ 1663 1611 ! cc $ sigmaw(i) 1664 1612 1665 1613 1666 ! ajout d'un effet onde de gravit é-Tgw(k)*deltatw(k) 03/02/06 YU1614 ! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU 1667 1615 ! Jingmei 1668 1616 … … 1676 1624 ! Sans GW 1677 1625 1678 ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)- spread(k)*deltatw(k))1626 ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-wkspread(k)*deltatw(k)) 1679 1627 1680 1628 ! GW formule 1 1681 1629 1682 1630 ! deltatw(k) = deltatw(k)+dtimesub* 1683 ! $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))1631 ! $ (ff+dtKE(k) - wkspread(k)*deltatw(k)-Tgw(k)*deltatw(k)) 1684 1632 1685 1633 ! GW formule 2 … … 1741 1689 ! Scale tendencies so that water vapour remains positive in w and x. 1742 1690 1743 CALL wake_vec_modulation(klon, klev, wk_adv, epsilon , qe, d_qe, deltaqw, &1691 CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qb, d_qb, deltaqw, & 1744 1692 d_deltaqw, sigmaw, d_sigmaw, alpha) 1693 ! 1694 ! Alpha_tot = Product of all the alpha's 1695 DO i = 1, klon 1696 IF (wk_adv(i)) THEN 1697 alpha_tot(i) = alpha_tot(i)*alpha(i) 1698 END IF 1699 END DO 1745 1700 1746 1701 ! cc nrlmd … … 1753 1708 DO i = 1, klon 1754 1709 IF (wk_adv(i) .AND. k<=kupper(i)) THEN 1755 d_t e(i, k) = alpha(i)*d_te(i, k)1756 d_q e(i, k) = alpha(i)*d_qe(i, k)1710 d_tb(i, k) = alpha(i)*d_tb(i, k) 1711 d_qb(i, k) = alpha(i)*d_qb(i, k) 1757 1712 d_deltatw(i, k) = alpha(i)*d_deltatw(i, k) 1758 1713 d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k) … … 1760 1715 #ifdef ISO 1761 1716 do ixt=1,ntraciso 1762 d_xt e(ixt,i,k)=alpha(i)*d_xte(ixt,i,k)1717 d_xtb(ixt,i,k)=alpha(i)*d_xtb(ixt,i,k) 1763 1718 d_deltaxtw(ixt,i,k)=alpha(i)*d_deltaxtw(ixt,i,k) 1764 1719 enddo !do ixt=1,ntraciso … … 1779 1734 DO i = 1, klon 1780 1735 IF (wk_adv(i) .AND. k<=kupper(i)) THEN 1781 dtls(i, k) = dtls(i, k) + d_t e(i, k)1782 dqls(i, k) = dqls(i, k) + d_q e(i, k)1736 dtls(i, k) = dtls(i, k) + d_tb(i, k) 1737 dqls(i, k) = dqls(i, k) + d_qb(i, k) 1783 1738 #ifdef ISO 1784 1739 do ixt=1,ntraciso 1785 dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xt e(ixt,i,k)1740 dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xtb(ixt,i,k) 1786 1741 enddo !do ixt=1,ntraciso 1787 1742 #endif … … 1817 1772 DO i = 1, klon 1818 1773 IF (wk_adv(i) .AND. k<=kupper(i)) THEN 1819 t e(i, k) = te0(i, k) + dtls(i, k)1820 q e(i, k) = qe0(i, k) + dqls(i, k)1821 th e(i, k) = te(i, k)/ppi(i, k)1774 tb(i, k) = tb0(i, k) + dtls(i, k) 1775 qb(i, k) = qb0(i, k) + dqls(i, k) 1776 thb(i, k) = tb(i, k)/ppi(i, k) 1822 1777 deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k) 1823 1778 deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k) 1824 1779 dth(i, k) = deltatw(i, k)/ppi(i, k) 1825 ! c print*,'k,qx,qw',k,q e(i,k)-sigmaw(i)*deltaqw(i,k)1826 ! c $ ,q e(i,k)+(1-sigmaw(i))*deltaqw(i,k)1780 ! c print*,'k,qx,qw',k,qb(i,k)-sigmaw(i)*deltaqw(i,k) 1781 ! c $ ,qb(i,k)+(1-sigmaw(i))*deltaqw(i,k) 1827 1782 #ifdef ISO 1828 1783 do ixt=1,ntraciso 1829 xt e(ixt,i,k) = xte0(ixt,i,k) + dxtls(ixt,i,k)1784 xtb(ixt,i,k) = xtb0(ixt,i,k) + dxtls(ixt,i,k) 1830 1785 deltaxtw(ixt,i,k) = deltaxtw(ixt,i,k)+d_deltaxtw(ixt,i,k) 1831 1786 enddo !do ixt=1,ntraciso … … 1848 1803 DO k= 1,klev 1849 1804 DO i = 1,klon 1850 call iso_verif_egalite_choix(xt e(iso_eau,i,k), &1851 q e(i,k),'wake 1379',errmax,errmaxrel)1805 call iso_verif_egalite_choix(xtb(iso_eau,i,k), & 1806 qb(i,k),'wake 1379',errmax,errmaxrel) 1852 1807 enddo ! DO i = 1,klon 1853 1808 enddo ! DO k= 1,klev … … 1855 1810 if (iso_hdo.gt.0) then 1856 1811 call iso_verif_aberrant_enc_vect2D( & 1857 xt e,qe, &1858 'wake 1456, xt eapres modifs',ntraciso,klon,klev)1812 xtb,qb, & 1813 'wake 1456, xtb apres modifs',ntraciso,klon,klev) 1859 1814 ! call iso_verif_aberrant_enc_vect2D_ns( 1860 1815 ! : deltaxtw,deltaqw, … … 1868 1823 DO i = 1, klon 1869 1824 IF (wk_adv(i)) THEN 1870 awdens(i) = awdens(i) + d_awdens(i) 1825 d_sig_gen2(i) = d_sig_gen2(i) + d_sig_gen(i) 1826 d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i) 1827 d_sig_col2(i) = d_sig_col2(i) + d_sig_col(i) 1828 d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i) 1829 d_sig_bnd2(i) = d_sig_bnd2(i) + d_sig_bnd(i) 1830 END IF 1831 END DO 1832 ! Bounds 1833 DO i = 1, klon 1834 IF (wk_adv(i)) THEN 1835 sigmaw_targ = max(sigmaw(i),sigmad) 1836 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) 1837 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 1838 sigmaw(i) = sigmaw_targ 1839 END IF 1840 END DO 1841 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1842 ! Cumulatives 1843 DO i = 1, klon 1844 IF (wk_adv(i)) THEN 1871 1845 wdens(i) = wdens(i) + d_wdens(i) 1872 d_awdens2(i) = d_awdens2(i) + d_awdens(i)1873 1846 d_wdens2(i) = d_wdens2(i) + d_wdens(i) 1847 d_dens_gen2(i) = d_dens_gen2(i) + d_dens_gen(i) 1848 d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i) 1849 d_dens_col2(i) = d_dens_col2(i) + d_dens_col(i) 1850 d_dens_bnd2(i) = d_dens_bnd2(i) + d_dens_bnd(i) 1874 1851 END IF 1875 1852 END DO 1853 ! Bounds 1876 1854 DO i = 1, klon 1877 1855 IF (wk_adv(i)) THEN 1878 1856 wdens_targ = max(wdens(i),wdensmin) 1857 d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i) 1879 1858 d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i) 1880 1859 wdens(i) = wdens_targ 1881 ! 1860 END IF 1861 END DO 1862 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1863 ! Cumulatives 1864 DO i = 1, klon 1865 IF (wk_adv(i)) THEN 1866 awdens(i) = awdens(i) + d_awdens(i) 1867 d_awdens2(i) = d_awdens2(i) + d_awdens(i) 1868 END IF 1869 END DO 1870 ! Bounds 1871 DO i = 1, klon 1872 IF (wk_adv(i)) THEN 1882 1873 wdens_targ = min( max(awdens(i),0.), wdens(i) ) 1874 d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) 1883 1875 d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i) 1884 1876 awdens(i) = wdens_targ 1885 1877 END IF 1886 1878 END DO 1887 DO i = 1, klon 1888 IF (wk_adv(i)) THEN 1889 sigmaw_targ = max(sigmaw(i),sigmad) 1890 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 1891 sigmaw(i) = sigmaw_targ 1892 END IF 1893 END DO 1879 ! 1880 IF (iflag_wk_pop_dyn >= 2) THEN 1881 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn >= 2!!!!!! 1882 ! Cumulatives 1883 DO i = 1, klon 1884 IF (wk_adv(i)) THEN 1885 d_adens_death2(i) = d_adens_death2(i) + d_adens_death(i) 1886 d_adens_icol2(i) = d_adens_icol2(i) + d_adens_icol(i) 1887 d_adens_acol2(i) = d_adens_acol2(i) + d_adens_acol(i) 1888 d_adens_bnd2(i) = d_adens_bnd2(i) + d_adens_bnd(i) 1889 END IF 1890 END DO 1891 ! Bounds 1892 DO i = 1, klon 1893 IF (wk_adv(i)) THEN 1894 wdens_targ = min( max(awdens(i),0.), wdens(i) ) 1895 d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) 1896 awdens(i) = wdens_targ 1897 END IF 1898 END DO 1899 ! 1900 IF (iflag_wk_pop_dyn == 3) THEN 1901 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! asigmaw for iflag_wk_pop_dyn = 3!!!!!! 1902 ! Cumulatives 1903 DO i = 1, klon 1904 IF (wk_adv(i)) THEN 1905 asigmaw(i) = asigmaw(i) + d_asigmaw(i) 1906 d_asigmaw2(i) = d_asigmaw2(i) + d_asigmaw(i) 1907 d_asig_death2(i) = d_asig_death2(i) + d_asig_death(i) 1908 d_asig_spread2(i) = d_asig_spread2(i) + d_asig_spread(i) 1909 d_asig_iicol2(i) = d_asig_iicol2(i) + d_asig_iicol(i) 1910 d_asig_aicol2(i) = d_asig_aicol2(i) + d_asig_aicol(i) 1911 d_asig_bnd2(i) = d_asig_bnd2(i) + d_asig_bnd(i) 1912 END IF 1913 END DO 1914 ! Bounds 1915 DO i = 1, klon 1916 IF (wk_adv(i)) THEN 1917 ! asigmaw lower bound set to sigmad/2 in order to allow asigmaw values lower than sigmad. 1918 !! sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i)) 1919 sigmaw_targ = min(max(asigmaw(i),sigmad/2.),sigmaw(i)) 1920 d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) 1921 d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) 1922 asigmaw(i) = sigmaw_targ 1923 END IF 1924 END DO 1925 1926 #ifdef IOPHYS_WK 1927 IF (phys_sub) THEN 1928 CALL iophys_ecrit('wdensb',1,'wdensb','m',wdens) 1929 CALL iophys_ecrit('awdensb',1,'awdensb','m',awdens) 1930 CALL iophys_ecrit('sigmawb',1,'sigmawb','m',sigmaw) 1931 CALL iophys_ecrit('asigmawb',1,'asigmawb','m',asigmaw) 1932 ! 1933 call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2) 1934 call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2) 1935 call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2) 1936 call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2) 1937 call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2) 1938 ! 1939 call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2) 1940 call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2) 1941 call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2) 1942 call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2) 1943 call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2) 1944 ! 1945 CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2) 1946 CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2) 1947 CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2) 1948 CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2) 1949 CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2) 1950 CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2) 1951 ! 1952 CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2) 1953 CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2) 1954 CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2) 1955 CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2) 1956 CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2) 1957 CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2) 1958 ENDIF 1959 #endif 1960 ENDIF ! (iflag_wk_pop_dyn == 3) 1961 ENDIF ! (iflag_wk_pop_dyn >= 2) 1894 1962 ENDIF ! (iflag_wk_pop_dyn >= 1) 1895 !>jyg 1896 1897 1898 ! Determine Ptop from buoyancy integral 1899 ! --------------------------------------- 1900 1901 ! - 1/ Pressure of the level where dth changes sign. 1902 1903 DO i = 1, klon 1904 IF (wk_adv(i)) THEN 1905 ptop_provis(i) = ph(i, 1) 1906 END IF 1907 END DO 1908 1909 DO k = 2, klev 1910 DO i = 1, klon 1911 IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & 1912 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 1913 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 1914 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 1915 END IF 1916 END DO 1917 END DO 1918 1919 ! - 2/ dth integral 1920 1921 DO i = 1, klon 1922 IF (wk_adv(i)) THEN !!! nrlmd 1923 sum_dth(i) = 0. 1924 dthmin(i) = -delta_t_min 1925 z(i) = 0. 1926 END IF 1927 END DO 1928 1929 DO k = 1, klev 1930 DO i = 1, klon 1931 IF (wk_adv(i)) THEN 1932 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) 1933 IF (dz(i)>0) THEN 1934 z(i) = z(i) + dz(i) 1935 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 1936 dthmin(i) = amin1(dthmin(i), dth(i,k)) 1937 END IF 1938 END IF 1939 END DO 1940 END DO 1941 1942 ! - 3/ height of triangle with area= sum_dth and base = dthmin 1943 1944 DO i = 1, klon 1945 IF (wk_adv(i)) THEN 1946 hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) 1947 hw(i) = amax1(hwmin, hw(i)) 1948 END IF 1949 END DO 1950 1951 ! - 4/ now, get Ptop 1952 1953 DO i = 1, klon 1954 IF (wk_adv(i)) THEN !!! nrlmd 1955 ktop(i) = 0 1956 z(i) = 0. 1957 END IF 1958 END DO 1959 1960 DO k = 1, klev 1961 DO i = 1, klon 1962 IF (wk_adv(i)) THEN 1963 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i)) 1964 IF (dz(i)>0) THEN 1965 z(i) = z(i) + dz(i) 1966 ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) 1967 ktop(i) = k 1968 END IF 1969 END IF 1970 END DO 1971 END DO 1972 1973 ! 4.5/Correct ktop and ptop 1974 1975 DO i = 1, klon 1976 IF (wk_adv(i)) THEN 1977 ptop_new(i) = ptop(i) 1978 END IF 1979 END DO 1980 1981 DO k = klev, 2, -1 1982 DO i = 1, klon 1983 ! IM v3JYG; IF (k .GE. ktop(i) 1984 IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & 1985 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 1986 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 1987 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 1988 END IF 1989 END DO 1990 END DO 1991 1992 1993 DO i = 1, klon 1994 IF (wk_adv(i)) THEN 1995 ptop(i) = ptop_new(i) 1996 END IF 1997 END DO 1998 1999 DO k = klev, 1, -1 2000 DO i = 1, klon 2001 IF (wk_adv(i)) THEN !!! nrlmd 2002 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 2003 END IF 2004 END DO 2005 END DO 1963 1964 1965 1966 Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, & 1967 dth, hw, rho, delta_t_min, & 1968 ktop, wk_adv, h_zzz, ptop1, ktop1) 1969 !! print'("pkupper APPEL ",7i6)',isubstep,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper 2006 1970 2007 1971 ! 5/ Set deltatw & deltaqw to 0 above kupper … … 2028 1992 DO i = 1, klon 2029 1993 IF (wk_adv(i)) THEN !!! nrlmd 2030 sum_th u(i) = 0.2031 sum_t u(i) = 0.2032 sum_q u(i) = 0.2033 sum_thv u(i) = 0.1994 sum_thx(i) = 0. 1995 sum_tx(i) = 0. 1996 sum_qx(i) = 0. 1997 sum_thvx(i) = 0. 2034 1998 sum_dth(i) = 0. 2035 1999 sum_dq(i) = 0. 2036 sum_rho(i) = 0.2037 2000 sum_dtdwn(i) = 0. 2038 2001 sum_dqdwn(i) = 0. 2039 2002 2040 av_th u(i) = 0.2041 av_t u(i) = 0.2042 av_q u(i) = 0.2043 av_thv u(i) = 0.2003 av_thx(i) = 0. 2004 av_tx(i) = 0. 2005 av_qx(i) = 0. 2006 av_thvx(i) = 0. 2044 2007 av_dth(i) = 0. 2045 2008 av_dq(i) = 0. 2046 av_rho(i) = 0.2047 2009 av_dtdwn(i) = 0. 2048 2010 av_dqdwn(i) = 0. … … 2053 2015 ! -------------------------------------- 2054 2016 2055 ! Initialize sum_thv uto 1st level virt. pot. temp.2017 ! Initialize sum_thvx to 1st level virt. pot. temp. 2056 2018 2057 2019 DO i = 1, klon … … 2059 2021 z(i) = 1. 2060 2022 dz(i) = 1. 2061 sum_thv u(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)2023 sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i) 2062 2024 sum_dth(i) = 0. 2063 2025 END IF … … 2067 2029 DO i = 1, klon 2068 2030 IF (wk_adv(i)) THEN !!! nrlmd 2069 dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)* rg)2031 dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG) 2070 2032 IF (dz(i)>0) THEN 2071 2033 z(i) = z(i) + dz(i) 2072 sum_th u(i) = sum_thu(i) + thu(i, k)*dz(i)2073 sum_t u(i) = sum_tu(i) + tu(i, k)*dz(i)2074 sum_q u(i) = sum_qu(i) + qu(i, k)*dz(i)2075 sum_thv u(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)2034 sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i) 2035 sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i) 2036 sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i) 2037 sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i) 2076 2038 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 2077 2039 sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) 2078 sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)2079 2040 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) 2080 2041 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) … … 2100 2061 DO i = 1, klon 2101 2062 IF (wk_adv(i)) THEN !!! nrlmd 2102 av_th u(i) = sum_thu(i)/hw0(i)2103 av_t u(i) = sum_tu(i)/hw0(i)2104 av_q u(i) = sum_qu(i)/hw0(i)2105 av_thv u(i) = sum_thvu(i)/hw0(i)2063 av_thx(i) = sum_thx(i)/hw0(i) 2064 av_tx(i) = sum_tx(i)/hw0(i) 2065 av_qx(i) = sum_qx(i)/hw0(i) 2066 av_thvx(i) = sum_thvx(i)/hw0(i) 2106 2067 av_dth(i) = sum_dth(i)/hw0(i) 2107 2068 av_dq(i) = sum_dq(i)/hw0(i) 2108 av_rho(i) = sum_rho(i)/hw0(i)2109 2069 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 2110 2070 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 2111 2071 2112 wape(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + & 2113 av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) 2114 END IF 2115 END DO 2072 wape(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + & 2073 av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i) 2074 END IF 2075 END DO 2076 2116 2077 2117 2078 ! Filter out bad wakes … … 2146 2107 !! sigmaw(i) = max(sigmad, sigd_con(i)) 2147 2108 sigmaw_targ = max(sigmad, sigd_con(i)) 2109 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) 2148 2110 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 2149 2111 sigmaw(i) = sigmaw_targ 2112 ! 2113 d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) 2114 d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) 2115 asigmaw(i) = sigmaw_targ 2150 2116 !>jyg 2151 2117 fip(i) = 0. … … 2157 2123 END IF 2158 2124 END DO 2159 2160 END DO ! end sub-timestep loop 2161 2125 ! 2126 ! ------------------------------------------------------------------------ 2127 ! 2128 END DO ! isubstep end sub-timestep loop 2129 ! 2130 ! ------------------------------------------------------------------------ 2131 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2132 ! ------------------------------------------------------------------------ 2133 ! 2134 2135 #ifdef IOPHYS_WK 2136 IF (.not.phys_sub) CALL iophys_ecrit('wape_b',1,'wape_b','J/kg',wape) 2137 #endif 2162 2138 IF (prt_level>=10) THEN 2163 2139 PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', & … … 2179 2155 ! cc 2180 2156 z(i) = 0. 2181 sum_th u(i) = 0.2182 sum_t u(i) = 0.2183 sum_q u(i) = 0.2184 sum_thv u(i) = 0.2157 sum_thx(i) = 0. 2158 sum_tx(i) = 0. 2159 sum_qx(i) = 0. 2160 sum_thvx(i) = 0. 2185 2161 sum_dth(i) = 0. 2186 2162 sum_half_dth(i) = 0. 2187 2163 sum_dq(i) = 0. 2188 sum_rho(i) = 0.2189 2164 sum_dtdwn(i) = 0. 2190 2165 sum_dqdwn(i) = 0. 2191 2166 2192 av_th u(i) = 0.2193 av_t u(i) = 0.2194 av_q u(i) = 0.2195 av_thv u(i) = 0.2167 av_thx(i) = 0. 2168 av_tx(i) = 0. 2169 av_qx(i) = 0. 2170 av_thvx(i) = 0. 2196 2171 av_dth(i) = 0. 2197 2172 av_dq(i) = 0. 2198 av_rho(i) = 0.2199 2173 av_dtdwn(i) = 0. 2200 2174 av_dqdwn(i) = 0. … … 2211 2185 IF (ok_qx_qw(i)) THEN 2212 2186 ! cc 2213 rho(i, k) = p(i, k)/( rd*te(i,k))2187 rho(i, k) = p(i, k)/(RD*tb(i,k)) 2214 2188 IF (k==1) THEN 2215 rhoh(i, k) = ph(i, k)/( rd*te(i,k))2189 rhoh(i, k) = ph(i, k)/(RD*tb(i,k)) 2216 2190 zhh(i, k) = 0 2217 2191 ELSE 2218 rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) 2219 zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) 2220 END IF 2221 the(i, k) = te(i, k)/ppi(i, k) 2222 thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 2223 tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) 2224 qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) 2225 rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) 2192 rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1))) 2193 zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1) 2194 END IF 2195 thb(i, k) = tb(i, k)/ppi(i, k) 2196 thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 2197 tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i) 2198 qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i) 2226 2199 dth(i, k) = deltatw(i, k)/ppi(i, k) 2227 2200 #ifdef ISO 2228 2201 do ixt=1,ntraciso 2229 xt u(ixt,i,k) = xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)2202 xtx(ixt,i,k) = xtb(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i) 2230 2203 enddo !do ixt=1,ntraciso 2231 2204 #endif … … 2238 2211 if (iso_hdo.gt.0) then 2239 2212 call iso_verif_aberrant_enc_vect2D( & 2240 xt u,qu, &2213 xtx,qx, & 2241 2214 'wake 1834, apres modifs',ntraciso,klon,klev) 2242 2215 endif … … 2247 2220 ! ----------------------------------------------------------- 2248 2221 2249 ! Initialize sum_thv uto 1st level virt. pot. temp.2222 ! Initialize sum_thvx to 1st level virt. pot. temp. 2250 2223 2251 2224 DO i = 1, klon … … 2256 2229 dz(i) = 1. 2257 2230 dz_half(i) = 1. 2258 sum_thv u(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)2231 sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i) 2259 2232 sum_dth(i) = 0. 2260 2233 END IF … … 2266 2239 IF (ok_qx_qw(i)) THEN 2267 2240 ! cc 2268 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)* rg)2269 dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)* rg)2241 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG) 2242 dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*RG) 2270 2243 IF (dz(i)>0) THEN 2271 2244 z(i) = z(i) + dz(i) 2272 sum_th u(i) = sum_thu(i) + thu(i, k)*dz(i)2273 sum_t u(i) = sum_tu(i) + tu(i, k)*dz(i)2274 sum_q u(i) = sum_qu(i) + qu(i, k)*dz(i)2275 sum_thv u(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)2245 sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i) 2246 sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i) 2247 sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i) 2248 sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i) 2276 2249 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 2277 2250 sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) 2278 sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)2279 2251 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) 2280 2252 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) … … 2306 2278 IF (ok_qx_qw(i)) THEN 2307 2279 ! cc 2308 av_th u(i) = sum_thu(i)/hw0(i)2309 av_t u(i) = sum_tu(i)/hw0(i)2310 av_q u(i) = sum_qu(i)/hw0(i)2311 av_thv u(i) = sum_thvu(i)/hw0(i)2280 av_thx(i) = sum_thx(i)/hw0(i) 2281 av_tx(i) = sum_tx(i)/hw0(i) 2282 av_qx(i) = sum_qx(i)/hw0(i) 2283 av_thvx(i) = sum_thvx(i)/hw0(i) 2312 2284 av_dth(i) = sum_dth(i)/hw0(i) 2313 2285 av_dq(i) = sum_dq(i)/hw0(i) 2314 av_rho(i) = sum_rho(i)/hw0(i)2315 2286 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 2316 2287 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 2317 2288 2318 wape2(i) = - rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &2319 av_dth(i)*av_q u(i)+av_dth(i)*av_dq(i)))/av_thvu(i)2289 wape2(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + & 2290 av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i) 2320 2291 END IF 2321 2292 END DO 2322 2293 #ifdef IOPHYS_WK 2294 IF (.not.phys_sub) CALL iophys_ecrit('wape2_a',1,'wape2_a','J/kg',wape2) 2295 #endif 2323 2296 2324 2297 … … 2350 2323 END DO 2351 2324 END IF 2325 #ifdef IOPHYS_WK 2326 IF (.not.phys_sub) CALL iophys_ecrit('wape2_b',1,'wape2_b','J/kg',wape2) 2327 #endif 2352 2328 2353 2329 … … 2384 2360 !! sigmaw(i) = amax1(sigmad, sigd_con(i)) 2385 2361 sigmaw_targ = max(sigmad, sigd_con(i)) 2362 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) 2386 2363 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 2387 2364 sigmaw(i) = sigmaw_targ 2365 ! 2366 d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) 2367 d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) 2368 asigmaw(i) = sigmaw_targ 2388 2369 !>jyg 2389 2370 fip(i) = 0. … … 2394 2375 gwake(i) = .TRUE. 2395 2376 END IF 2396 END IF 2377 #ifdef IOPHYS_WK 2378 IF (.not.phys_sub) CALL iophys_ecrit('cstar2',1,'cstar2','J/kg',cstar2) 2379 #endif 2380 END IF ! (ok_qx_qw(i)) 2397 2381 END DO 2398 2382 … … 2426 2410 END IF 2427 2411 END DO 2428 2412 IF (iflag_wk_pop_dyn >= 3) THEN 2413 #ifdef IOPHYS_WK 2414 IF (.not.phys_sub) THEN 2415 CALL iophys_ecrit('fip',1,'fip','J/kg',fip) 2416 CALL iophys_ecrit('hw',1,'hw','J/kg',hw) 2417 CALL iophys_ecrit('ptop',1,'ptop','J/kg',ptop) 2418 CALL iophys_ecrit('wdens',1,'wdens','J/kg',wdens) 2419 CALL iophys_ecrit('awdens',1,'awdens','m',awdens) 2420 CALL iophys_ecrit('sigmaw',1,'sigmaw','m',sigmaw) 2421 CALL iophys_ecrit('asigmaw',1,'asigmaw','m',asigmaw) 2422 ! 2423 CALL iophys_ecrit('rad_wk',1,'rad_wk','J/kg',rad_wk) 2424 CALL iophys_ecrit('arad_wk',1,'arad_wk','J/kg',arad_wk) 2425 CALL iophys_ecrit('irad_wk',1,'irad_wk','J/kg',irad_wk) 2426 ! 2427 call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2) 2428 call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2) 2429 call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2) 2430 call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2) 2431 call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2) 2432 ! 2433 call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2) 2434 call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2) 2435 call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2) 2436 call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2) 2437 call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2) 2438 ! 2439 CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2) 2440 CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2) 2441 CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2) 2442 CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2) 2443 CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2) 2444 CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2) 2445 ! 2446 CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2) 2447 CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2) 2448 CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2) 2449 CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2) 2450 CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2) 2451 CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2) 2452 ENDIF ! (.not.phys_sub) 2453 #endif 2454 ENDIF ! (iflag_wk_pop_dyn >= 3) 2429 2455 ! Limitation de sigmaw 2430 2456 … … 2441 2467 DO i = 1, klon 2442 2468 kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & 2443 .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin) 2469 .NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold) 2470 !! .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin) 2444 2471 ENDDO 2445 2472 ELSE ! (iflag_wk_pop_dyn >= 1) … … 2481 2508 wape(i) = 0. 2482 2509 cstar(i) = 0. 2483 !!jyg Outside subroutine "Wake" hw, wdens andsigmaw are zero when there are no wakes2510 !!jyg Outside subroutine "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes 2484 2511 !! hw(i) = hwmin !jyg 2485 2512 !! sigmaw(i) = sigmad !jyg 2486 2513 hw(i) = 0. !jyg 2487 2514 fip(i) = 0. 2515 ! 2488 2516 !! sigmaw(i) = 0. !jyg 2489 2517 sigmaw_targ = 0. 2490 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 2518 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) 2519 !! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 2520 d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i) ! _in = correction jyg 20220124 2491 2521 sigmaw(i) = sigmaw_targ 2522 ! 2523 IF (iflag_wk_pop_dyn >= 3) THEN 2524 sigmaw_targ = 0. 2525 d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) 2526 !! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 2527 d_asigmaw2(i) = sigmaw_targ - asigmaw_in(i) ! _in = correction jyg 20220124 2528 asigmaw(i) = sigmaw_targ 2529 ELSE 2530 asigmaw(i) = 0. 2531 ENDIF ! (iflag_wk_pop_dyn >= 3) 2532 ! 2492 2533 IF (iflag_wk_pop_dyn >= 1) THEN 2493 2534 !! awdens(i) = 0. 2494 2535 !! wdens(i) = 0. 2495 2536 wdens_targ = 0. 2496 d_wdens2(i) = wdens_targ - wdens(i) 2537 d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i) 2538 !! d_wdens2(i) = wdens_targ - wdens(i) 2539 d_wdens2(i) = wdens_targ - wdens_in(i) ! jyg 20220916 2497 2540 wdens(i) = wdens_targ 2498 2541 wdens_targ = 0. 2499 d_awdens2(i) = wdens_targ - awdens(i) 2542 !!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens. 2543 IF (iflag_wk_pop_dyn >= 2) THEN 2544 d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) 2545 ENDIF ! (iflag_wk_pop_dyn >= 2) 2546 !! d_awdens2(i) = wdens_targ - awdens(i) 2547 d_awdens2(i) = wdens_targ - awdens_in(i) ! jyg 20220916 2500 2548 awdens(i) = wdens_targ 2549 !! IF (iflag_wk_pop_dyn == 2) THEN 2550 !! d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) 2551 !! ENDIF ! (iflag_wk_pop_dyn == 2) 2501 2552 ENDIF ! (iflag_wk_pop_dyn >= 1) 2502 2553 ELSE ! (kill_wake(i)) … … 2512 2563 wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout) 2513 2564 ENDIF 2565 #ifdef IOPHYS_WK 2566 IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape) 2567 #endif 2514 2568 2515 2569 … … 2543 2597 END DO 2544 2598 !jyg< 2545 DO i = 1, klon 2546 d_sigmaw2(i) = d_sigmaw2(i)/dtime 2547 d_awdens2(i) = d_awdens2(i)/dtime 2548 d_wdens2(i) = d_wdens2(i)/dtime 2549 ENDDO 2599 IF (iflag_wk_pop_dyn >= 1) THEN 2600 DO i = 1, klon 2601 IF (ok_qx_qw(i)) THEN 2602 d_sig_gen2(i) = d_sig_gen2(i)/dtime 2603 d_sig_death2(i) = d_sig_death2(i)/dtime 2604 d_sig_col2(i) = d_sig_col2(i)/dtime 2605 d_sig_spread2(i) = d_sig_spread2(i)/dtime 2606 d_sig_bnd2(i) = d_sig_bnd2(i)/dtime 2607 d_sigmaw2(i) = d_sigmaw2(i)/dtime 2608 ! 2609 d_dens_gen2(i) = d_dens_gen2(i)/dtime 2610 d_dens_death2(i) = d_dens_death2(i)/dtime 2611 d_dens_col2(i) = d_dens_col2(i)/dtime 2612 d_dens_bnd2(i) = d_dens_bnd2(i)/dtime 2613 d_awdens2(i) = d_awdens2(i)/dtime 2614 d_wdens2(i) = d_wdens2(i)/dtime 2615 ENDIF 2616 ENDDO 2617 IF (iflag_wk_pop_dyn >= 2) THEN 2618 DO i = 1, klon 2619 IF (ok_qx_qw(i)) THEN 2620 d_adens_death2(i) = d_adens_death2(i)/dtime 2621 d_adens_icol2(i) = d_adens_icol2(i)/dtime 2622 d_adens_acol2(i) = d_adens_acol2(i)/dtime 2623 d_adens_bnd2(i) = d_adens_bnd2(i)/dtime 2624 ENDIF 2625 ENDDO 2626 IF (iflag_wk_pop_dyn == 3) THEN 2627 DO i = 1, klon 2628 IF (ok_qx_qw(i)) THEN 2629 d_asig_death2(i) = d_asig_death2(i)/dtime 2630 d_asig_iicol2(i) = d_asig_iicol2(i)/dtime 2631 d_asig_aicol2(i) = d_asig_aicol2(i)/dtime 2632 d_asig_spread2(i) = d_asig_spread2(i)/dtime 2633 d_asig_bnd2(i) = d_asig_bnd2(i)/dtime 2634 ENDIF 2635 ENDDO 2636 ENDIF ! (iflag_wk_pop_dyn == 3) 2637 ENDIF ! (iflag_wk_pop_dyn >= 2) 2638 ENDIF ! (iflag_wk_pop_dyn >= 1) 2639 2550 2640 !>jyg 2551 2641 2552 2553 2554 RETURN 2642 RETURN 2555 2643 END SUBROUTINE wake 2556 2644 2557 SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon , qe, d_qe, deltaqw, &2645 SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qb, d_qb, deltaqw, & 2558 2646 d_deltaqw, sigmaw, d_sigmaw, alpha) 2559 2647 ! ------------------------------------------------------ 2560 ! D \'etermination du coefficient alpha tel que les tendances2648 ! Dtermination du coefficient alpha tel que les tendances 2561 2649 ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent 2562 2650 ! a une humidite positive dans la zone (x) et dans la zone (w). … … 2565 2653 2566 2654 ! Input 2567 REAL q e(nlon, nl), d_qe(nlon, nl)2655 REAL qb(nlon, nl), d_qb(nlon, nl) 2568 2656 REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl) 2569 2657 REAL sigmaw(nlon), d_sigmaw(nlon) … … 2576 2664 REAL alpha1(nlon) 2577 2665 REAL x, a, b, c, discrim 2578 REAL epsilon 2579 ! DATA epsilon/1.e-15/ 2666 REAL epsilon_loc 2580 2667 INTEGER i,k 2581 2668 … … 2592 2679 DO i = 1, nlon 2593 2680 IF (wk_adv(i)) THEN 2594 x = q e(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &2681 x = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qb(i, k) + & 2595 2682 (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * & 2596 2683 (deltaqw(i,k)+d_deltaqw(i,k)) 2597 2684 a = -d_sigmaw(i)*d_deltaqw(i, k) 2598 b = d_q e(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &2685 b = d_qb(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - & 2599 2686 deltaqw(i, k)*d_sigmaw(i) 2600 c = q e(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon2687 c = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc 2601 2688 discrim = b*b - 4.*a*c 2602 2689 ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim 2603 IF (a+b>=0.) THEN !! Condition suffisante pour la positivit éde ovap2690 IF (a+b>=0.) THEN !! Condition suffisante pour la positivite de ovap 2604 2691 alpha1(i) = 1. 2605 2692 ELSE … … 2627 2714 END SUBROUTINE wake_vec_modulation 2628 2715 2716 2717 2718 SUBROUTINE pkupper (klon, klev, ptop, ph, p, pupper, kupper, & 2719 dth, hw_, rho, delta_t_min_in, & 2720 ktop, wk_adv, h_zzz, ptop1, ktop1) 2721 2722 USE lmdz_wake_ini , ONLY : wk_pupper 2723 USE lmdz_wake_ini , ONLY : RG 2724 USE lmdz_wake_ini , ONLY : hwmin 2725 USE lmdz_wake_ini , ONLY : iflag_wk_new_ptop, wk_delta_t_min, wk_frac_int_delta_t 2726 USE lmdz_wake_ini , ONLY : wk_int_delta_t_min 2727 2728 IMPLICIT NONE 2729 2730 INTEGER, INTENT(IN) :: klon,klev 2731 REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: ph, p 2732 REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: rho 2733 LOGICAL, DIMENSION (klon) , INTENT(IN) :: wk_adv 2734 REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: dth 2735 REAL, INTENT(IN) :: delta_t_min_in 2736 2737 2738 REAL, DIMENSION (klon) , INTENT(OUT) :: hw_ 2739 REAL, DIMENSION (klon) , INTENT(OUT) :: ptop 2740 INTEGER, DIMENSION (klon) , INTENT(OUT) :: Ktop 2741 REAL, DIMENSION (klon) , INTENT(OUT) :: pupper 2742 INTEGER, DIMENSION (klon) , INTENT(OUT) :: kupper 2743 REAL, DIMENSION (klon) , INTENT(OUT) :: h_zzz !! 2744 REAL, DIMENSION (klon) , INTENT(OUT) :: Ptop1 !! 2745 INTEGER, DIMENSION (klon) , INTENT(OUT) :: ktop1 !! 2746 2747 INTEGER :: i,k 2748 2749 LOGICAL, DIMENSION (klon) :: wk_active 2750 REAL :: delta_t_min 2751 REAL, DIMENSION (klon) :: dthmin 2752 REAL, DIMENSION (klon) :: ptop_provis,ptop_new 2753 REAL, DIMENSION (klon) :: z, dz 2754 REAL, DIMENSION (klon) :: sum_dth 2755 2756 INTEGER, DIMENSION (klon) :: k_ptop_provis 2757 REAL, DIMENSION (klon) :: zk_ptop_provis 2758 REAL, DIMENSION (klon) :: omega !! 2759 REAL, DIMENSION (klon,klev+1) :: int_dth !! 2760 REAL, DIMENSION (klon,klev+1) :: dzz !! 2761 REAL, DIMENSION (klon,klev+1) :: zzz !! 2762 REAL, DIMENSION (klon) :: frac_int_dth !! 2763 REAL :: ddd!! 2764 2765 2766 INTEGER, SAVE :: ipas=0 2767 2768 2769 2770 !INTEGER, SAVE :: compte=0 2771 2772 ! LJYF : a priori z, dz sum_dth sont aussi des variables internes 2773 ! Les eliminer apres verification convergence numerique 2774 2775 !compte=compte+1 2776 !print*,'compte=',compte 2777 2778 ! Determine Ptop from buoyancy integral 2779 ! --------------------------------------- 2780 2781 ! - 1/ Pressure of the level where dth changes sign. 2782 !print*,'WAKE LJYF' 2783 2784 2785 if (iflag_wk_new_ptop==0) then 2786 delta_t_min=delta_t_min_in 2787 else 2788 delta_t_min=wk_delta_t_min 2789 endif 2790 2791 DO i = 1, klon 2792 ptop_provis(i) = ph(i, 1) 2793 k_ptop_provis(i) = 1 2794 END DO 2795 2796 DO k = 2, klev 2797 DO i = 1, klon 2798 IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & 2799 ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2800 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2801 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 2802 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 2803 k_ptop_provis(i) = k 2804 END IF 2805 END DO 2806 END DO 2807 2808 2809 2810 ! - 2/ dth integral 2811 2812 DO i = 1, klon 2813 IF (wk_adv(i)) THEN !!! nrlmd 2814 sum_dth(i) = 0. 2815 dthmin(i) = -delta_t_min 2816 z(i) = 0. 2817 END IF 2818 END DO 2819 2820 DO k = 1, klev 2821 DO i = 1, klon 2822 IF (wk_adv(i)) THEN 2823 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG) 2824 IF (dz(i)>0) THEN 2825 z(i) = z(i) + dz(i) 2826 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 2827 dthmin(i) = amin1(dthmin(i), dth(i,k)) 2828 END IF 2829 END IF 2830 END DO 2831 END DO 2832 2833 ! - 3/ height of triangle with area= sum_dth and base = dthmin 2834 2835 DO i = 1, klon 2836 IF (wk_adv(i)) THEN 2837 hw_(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) 2838 hw_(i) = amax1(hwmin, hw_(i)) 2839 END IF 2840 END DO 2841 2842 ! - 4/ now, get Ptop 2843 2844 DO i = 1, klon 2845 IF (wk_adv(i)) THEN !!! nrlmd 2846 ktop(i) = 0 2847 z(i) = 0. 2848 END IF 2849 END DO 2850 2851 DO k = 1, klev 2852 DO i = 1, klon 2853 IF (wk_adv(i)) THEN 2854 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw_(i)-z(i)) 2855 IF (dz(i)>0) THEN 2856 z(i) = z(i) + dz(i) 2857 ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i) 2858 ktop(i) = k 2859 END IF 2860 END IF 2861 END DO 2862 END DO 2863 2864 ! 4.5/Correct ktop and ptop 2865 2866 DO i = 1, klon 2867 ptop_new(i) = ptop(i) 2868 END DO 2869 2870 DO k = klev, 2, -1 2871 DO i = 1, klon 2872 ! IM v3JYG; IF (k .GE. ktop(i) 2873 IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & 2874 ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2875 dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2876 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 2877 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 2878 END IF 2879 END DO 2880 END DO 2881 2882 2883 DO i = 1, klon 2884 ptop(i) = ptop_new(i) 2885 END DO 2886 2887 DO k = klev, 1, -1 2888 DO i = 1, klon 2889 IF (wk_adv(i)) THEN !!! nrlmd 2890 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 2891 END IF 2892 END DO 2893 END DO 2894 2895 ! IF (prt_level>=10) THEN 2896 ! PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) 2897 ! ENDIF 2898 2899 ! ----------------------------------------------------------------------- 2900 ! nouveau calcul de hw et ptop 2901 ! ----------------------------------------------------------------------- 2902 !if (iflag_wk_new_ptop>0) then 2903 do i=1,klon 2904 ptop1(i)=ph(i,1) 2905 ktop1(i)=1 2906 h_zzz(i)=0. 2907 enddo 2908 2909 IF (iflag_wk_new_ptop/=0) THEN 2910 2911 int_dth(1:klon,1:klev+1)=0. 2912 DO i = 1, klon 2913 IF (wk_adv(i)) THEN 2914 int_dth(i,1) = 0. 2915 END IF 2916 END DO 2917 2918 if (abs(iflag_wk_new_ptop) == 1 ) then 2919 DO k = 2, klev+1 2920 Do i = 1, klon 2921 IF (wk_adv(i)) THEN 2922 if (k<=k_ptop_provis(i)) then 2923 ddd=dth(i,k-1)*(ph(i,k-1) - max(ptop_provis(i),ph(i,k))) 2924 !ddd=dth(i,k-1)*(ph(i,k-1) - ph(i,k)) 2925 else 2926 ddd=0. 2927 endif 2928 int_dth(i,k) = int_dth(i,k-1) + ddd 2929 !ELSE 2930 ! int_dth(i,k) = 0. 2931 END IF 2932 END DO 2933 END DO 2934 else 2935 k_ptop_provis(:)=klev+1 2936 dthmin(:)=dth(:,1) 2937 ! calcul de l'int??grale de dT * dP jusqu'au dernier 2938 ! niveau avec dT<0. (en s'assurant qu'on a bien un 2939 ! dT negatif plus bas) 2940 DO k = 1, klev 2941 DO i = 1, klon 2942 dthmin(i)=min(dthmin(i),dth(i,k)) 2943 ddd=dth(i,k)*(ph(i,k)-ph(i,k+1)) 2944 if (dthmin(i)<0.) then 2945 if (k>=k_ptop_provis(i)) then 2946 ddd=0. 2947 else if (dth(i,k)>=0.) then 2948 ddd=0. 2949 k_ptop_provis(i)=k+1 2950 endif 2951 endif 2952 int_dth(i,k+1) = int_dth(i,k)+ ddd 2953 ENDDO 2954 ENDDO 2955 2956 DO i = 1, klon 2957 if ( k_ptop_provis(i)==klev+1 .or. .not. wk_adv(i)) then 2958 k_ptop_provis(i)=1 2959 endif 2960 ENDDO 2961 endif ! (abs(iflag_wk_new_ptop) == 1 ) 2962 ! print*, 'xxx, int_dth', (k,int_dth(1,k),k=1,klev) 2963 ! print*, 'xxx, k_ptop_provis', k_ptop_provis(1) 2964 2965 2966 2967 ! On se limite ?? des poches avec integrale dT * dp < -wk_int_delta_t_min 2968 do i=1,klon 2969 if (int_dth(i,k_ptop_provis(i)) > -wk_int_delta_t_min .or. k_ptop_provis(i)==1) then 2970 !if (1==0) then 2971 wk_active(i)=.false. 2972 ptop(i)=ph(i,1) 2973 ktop(i)=1 2974 hw_(i)=0. 2975 else 2976 wk_active(i)=wk_adv(i) 2977 endif 2978 enddo 2979 2980 DO i=1,klon 2981 IF (wk_active(i)) THEN 2982 frac_int_dth(i)=wk_frac_int_delta_t*int_dth(i,k_ptop_provis(i)) 2983 ENDIF 2984 ENDDO 2985 DO k = 1,klev 2986 DO i =1, klon 2987 ! print*,ipas,'yyy ',k,int_dth(i,k),frac_int_dth(i) 2988 IF (wk_active(i)) THEN 2989 IF (int_dth(i,k)>=frac_int_dth(i)) THEN 2990 ktop1(i) = min(k, k_ptop_provis(i)) 2991 !ktop1(i) = k 2992 !print*,ipas,'yyy ktop1= ',ktop1 2993 ENDIF 2994 ENDIF 2995 END DO 2996 END DO 2997 !print*, 'LAMINE' 2998 2999 DO i = 1, klon 3000 IF (wk_active(i)) THEN 3001 !print*, ipas,'xxx1, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ',ktop1 3002 ddd=int_dth(i,ktop1(i)+1)-int_dth(i,ktop1(i)) 3003 if (ddd==0.) then 3004 omega(i)=0. 3005 else 3006 omega(i) = (frac_int_dth(i) - int_dth(i,ktop1(i)))/ddd 3007 endif 3008 !! print*,'OMEGA ',omega(i) 3009 END IF 3010 END DO 3011 3012 !! print*, 'xxx' 3013 DO i = 1, klon 3014 IF (wk_active(i)) THEN 3015 ! print*, 'xxx, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ', & 3016 ! int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) 3017 ! print*, 'xxx, omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ', & 3018 !e omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) 3019 ptop1(i) = min((1 - omega(i))*ph(i,ktop1(i)) + omega(i)*ph(i,ktop1(i)+1), ph(i,1)) 3020 END IF 3021 END DO 3022 3023 DO i=1, klon 3024 IF (wk_active(i)) THEN 3025 zzz(i, 1) = 0 3026 END IF 3027 END DO 3028 DO k = 1, klev 3029 DO i = 1, klon 3030 IF (wk_active(i)) THEN 3031 dzz(i,k) = (ph(i,k) - ph(i,k+1))/(rho(i,k)*RG) 3032 zzz(i,k+1) = zzz(i,k) + dzz(i,k) 3033 END IF 3034 END DO 3035 END DO 3036 3037 DO i =1, klon 3038 IF (wk_active(i)) THEN 3039 h_zzz(i) = max((1- omega(i))*zzz(i,ktop1(i)) + omega(i)*zzz(i,ktop1(i)+1), hwmin) 3040 END IF 3041 END DO 3042 3043 3044 ENDIF ! (iflag_wk_new_ptop/=0) 3045 3046 !if (iflag_wk_new_ptop==2) then 3047 IF (iflag_wk_new_ptop>0) THEN 3048 do i=1,klon 3049 ptop(i)=ptop1(i) 3050 ktop(i)=ktop1(i) 3051 hw_(i)=h_zzz(i) 3052 enddo 3053 3054 !endif 3055 ENDIF 3056 3057 kupper = 0 3058 3059 IF (wk_pupper<1.) THEN 3060 ! Choose an integration bound well above wake top 3061 ! ----------------------------------------------------------------- 3062 3063 ! Pupper = 50000. ! melting level 3064 ! Pupper = 60000. 3065 ! Pupper = 80000. ! essais pour case_e 3066 DO i = 1, klon 3067 ! pupper(i) = 0.6*ph(i, 1) 3068 pupper(i) = wk_pupper*ph(i, 1) 3069 pupper(i) = max(pupper(i), 45000.) 3070 ! cc Pupper(i) = 60000. 3071 END DO 3072 3073 ELSE 3074 DO i=1, klon 3075 ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) 3076 ! pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-50.) 3077 pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-5000.) 3078 END DO 3079 END IF 3080 3081 ! -5/ Determination de kupper 3082 3083 DO k = klev, 1, -1 3084 DO i = 1, klon 3085 IF (ph(i,k+1)<pupper(i)) kupper(i) = k 3086 END DO 3087 END DO 3088 3089 ! On evite kupper = 1 et kupper = klev 3090 DO i = 1, klon 3091 kupper(i) = max(kupper(i), 2) 3092 kupper(i) = min(kupper(i), klev-1) 3093 END DO 3094 !---------- FIN nouveau calcul hw et ptop ------------------------------------- 3095 3096 IF (iflag_wk_new_ptop==999) THEN 3097 DO i = 1, klon 3098 hw_(i)=0. 3099 ptop(i)=ph(i,1) 3100 Ktop(i)=1 3101 pupper(i)=ph(i,2) 3102 kupper(i)=2 3103 h_zzz(i)=0. 3104 Ptop1(i)=ph(i,1) 3105 ENDDO 3106 ENDIF 3107 3108 zk_ptop_provis=k_ptop_provis 3109 3110 RETURN 3111 END SUBROUTINE pkupper 3112 3113 3114 SUBROUTINE wake_popdyn_1(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, & 3115 wdensmin, & 3116 dtimesub, gfl, rad_wk, f_shear, drdt_pos, & 3117 d_awdens, d_wdens, d_sigmaw, & 3118 iflag_wk_act, wk_adv, cin, wape, & 3119 drdt, & 3120 d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & 3121 d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & 3122 d_wdens_targ, d_sigmaw_targ) 3123 3124 3125 USE lmdz_wake_ini , ONLY : wake_ini 3126 USE lmdz_wake_ini , ONLY : prt_level,RG 3127 USE lmdz_wake_ini , ONLY : stark, wdens_ref 3128 USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0 3129 !! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin 3130 USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn 3131 USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max 3132 3133 IMPLICIT NONE 3134 3135 INTEGER, INTENT(IN) :: klon,klev 3136 LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv 3137 REAL, INTENT(IN) :: dtime 3138 REAL, INTENT(IN) :: dtimesub 3139 REAL, INTENT(IN) :: wdensmin 3140 REAL, DIMENSION (klon), INTENT(IN) :: wgen 3141 REAL, DIMENSION (klon), INTENT(IN) :: wdens 3142 REAL, DIMENSION (klon), INTENT(IN) :: awdens 3143 REAL, DIMENSION (klon), INTENT(IN) :: sigmaw 3144 REAL, DIMENSION (klon), INTENT(IN) :: cstar 3145 REAL, DIMENSION (klon), INTENT(IN) :: cin, wape 3146 REAL, DIMENSION (klon), INTENT(IN) :: f_shear 3147 INTEGER, INTENT(IN) :: iflag_wk_act 3148 3149 3150 ! 3151 3152 ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields 3153 ! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens) 3154 REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk 3155 REAL, DIMENSION (klon), INTENT(OUT) :: gfl 3156 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_awdens, d_wdens 3157 REAL, DIMENSION (klon), INTENT(OUT) :: drdt 3158 ! Some components of the tendencies of state variables 3159 REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd 3160 REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_spread 3161 REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd 3162 REAL, INTENT(OUT) :: d_wdens_targ, d_sigmaw_targ 3163 3164 3165 REAL :: delta_t_min 3166 INTEGER :: i, k 3167 REAL :: wdens0 3168 ! IM 080208 3169 LOGICAL, DIMENSION (klon) :: gwake 3170 3171 ! Variables liees a la dynamique de population 3172 REAL, DIMENSION(klon) :: act 3173 REAL, DIMENSION(klon) :: tau_wk_inv 3174 REAL, DIMENSION(klon) :: wape1_act, wape2_act 3175 LOGICAL, DIMENSION (klon) :: kill_wake 3176 REAL :: drdt_pos 3177 REAL :: tau_wk_inv_min 3178 3179 3180 3181 IF (iflag_wk_act == 0) THEN 3182 act(:) = 0. 3183 ELSEIF (iflag_wk_act == 1) THEN 3184 act(:) = 1. 3185 ELSEIF (iflag_wk_act ==2) THEN 3186 DO i = 1, klon 3187 IF (wk_adv(i)) THEN 3188 wape1_act(i) = abs(cin(i)) 3189 wape2_act(i) = 2.*wape1_act(i) + 1. 3190 act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) )) 3191 ENDIF ! (wk_adv(i)) 3192 ENDDO 3193 ENDIF ! (iflag_wk_act ==2) 3194 3195 DO i = 1, klon 3196 IF (wk_adv(i)) THEN 3197 rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero) 3198 gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) 3199 END IF 3200 END DO 3201 3202 DO i = 1, klon 3203 IF (wk_adv(i)) THEN 3204 !! tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.) 3205 tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.) 3206 tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub) 3207 drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / & 3208 (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i)) 3209 !! (1 - 2*sigmaw(i)*(1.-f_shear(i))) 3210 drdt_pos=max(drdt(i),0.) 3211 3212 !! d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) & 3213 !! - wdens(i)*tau_wk_inv_min & 3214 !! - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub 3215 !jyg+mlt< 3216 d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub 3217 d_dens_gen(i) = wgen(i) 3218 d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min 3219 d_dens_col(i) = -2.*wdens(i)*gfl(i)*drdt_pos 3220 d_dens_gen(i) = d_dens_gen(i)*dtimesub 3221 d_dens_death(i) = d_dens_death(i)*dtimesub 3222 d_dens_col(i) = d_dens_col(i)*dtimesub 3223 3224 d_wdens(i) = d_dens_gen(i)+d_dens_death(i)+d_dens_col(i) 3225 !! d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min - & 3226 !! 2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub 3227 !>jyg+mlt 3228 ! 3229 !jyg< 3230 d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i)) 3231 !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) 3232 d_dens_bnd(i) = d_wdens_targ - d_wdens(i) 3233 d_wdens(i) = d_wdens_targ 3234 !! d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i)) 3235 !>jyg 3236 3237 !jyg+mlt< 3238 !! d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) & 3239 !! + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) & 3240 !! - sigmaw(i)*tau_wk_inv_min )*dtimesub 3241 d_sig_gen(i) = wgen(i)*aa0 3242 d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min 3243 !! 3244 3245 d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos 3246 d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos 3247 d_sig_spread(i) = gfl(i)*cstar(i) 3248 d_sig_gen(i) = d_sig_gen(i)*dtimesub 3249 d_sig_death(i) = d_sig_death(i)*dtimesub 3250 d_sig_col(i) = d_sig_col(i)*dtimesub 3251 d_sig_spread(i) = d_sig_spread(i)*dtimesub 3252 d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) 3253 !>jyg+mlt 3254 ! 3255 !jyg< 3256 d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i)) 3257 !! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) 3258 !! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) 3259 d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) 3260 d_sigmaw(i) = d_sigmaw_targ 3261 !! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) 3262 !>jyg 3263 ENDIF 3264 ENDDO 3265 3266 IF (prt_level >= 10) THEN 3267 print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', & 3268 cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) 3269 print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', & 3270 wdens(1), awdens(1), act(1), d_awdens(1) 3271 print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', & 3272 wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1) 3273 print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & 3274 d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) 3275 ENDIF 3276 3277 RETURN 3278 END SUBROUTINE wake_popdyn_1 3279 3280 SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, & 3281 wdensmin, & 3282 sigmaw, wdens, awdens, & !! states variables 3283 gfl, cstar, cin, wape, rad_wk, & 3284 d_sigmaw, d_wdens, d_awdens, & !! tendences 3285 cont_fact, & 3286 d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & 3287 d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & 3288 d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd ) 3289 3290 3291 3292 USE lmdz_wake_ini , ONLY : wake_ini 3293 USE lmdz_wake_ini , ONLY : prt_level,RG 3294 USE lmdz_wake_ini , ONLY : stark, wdens_ref 3295 USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0 3296 !! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin 3297 USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn 3298 USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max 3299 3300 IMPLICIT NONE 3301 3302 INTEGER, INTENT(IN) :: klon,klev 3303 LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv 3304 REAL, INTENT(IN) :: dtimesub 3305 REAL, INTENT(IN) :: wdensmin 3306 REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes 3307 REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes 3308 REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area 3309 REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area 3310 REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes 3311 REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre 3312 3313 ! 3314 REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = wake radius 3315 REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front lenght per unit area 3316 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_wdens, d_awdens 3317 REAL, DIMENSION (klon), INTENT(OUT) :: cont_fact !! RM facteur de contact = 2 pi * rad * C* 3318 ! Some components of the tendencies of state variables 3319 REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd 3320 REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd 3321 REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd 3322 3323 3324 !! internal variables 3325 3326 INTEGER :: i, k 3327 REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes 3328 REAL :: tau_wk_inv_min 3329 REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes 3330 REAL :: d_wdens_targ, d_sigmaw_targ 3331 3332 3333 !! Equations 3334 !! dD/dt = B - (D-A)/tau - f D^2 3335 !! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^2 3336 !! dsigma/dt = B a0 - sigma/D (D-A)/tau + Lg C* - f (D-A)^2 (sigma/D-a0) 3337 !! 3338 !! f = 2 (B (a0-sigma/D) + Lg C*) / (2 (D-A)^2 (2 sigma/D-a0) + D (1-2 sigma)) 3339 3340 3341 DO i = 1, klon 3342 IF (wk_adv(i)) THEN 3343 rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero) 3344 gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) 3345 END IF 3346 END DO 3347 3348 3349 DO i = 1, klon 3350 IF (wk_adv(i)) THEN 3351 !! tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.) 3352 tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.) 3353 tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub) 3354 tau_prime(i) = tau_cv 3355 !! cont_fact(i) = 2.*(wgen(i)*(aa0-sigmaw(i)/wdens(i)) + gfl(i)*cstar(i)) / & 3356 !! (2.*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i) - aa0) + wdens(i)*(1.-2.*sigmaw(i))) 3357 !! cont_fact(i) = 2.*3.14*rad_wk(i)*cstar(i) ! bug 3358 !! cont_fact(i) = 4.*3.14*rad_wk(i)*cstar(i) 3359 cont_fact(i) = 2.*gfl(i)*cstar(i)/wdens(i) 3360 3361 d_sig_gen(i) = wgen(i)*aa0 3362 d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min 3363 d_sig_col(i) = - cont_fact(i)*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i)-aa0) 3364 d_sig_spread(i) = gfl(i)*cstar(i) 3365 ! 3366 d_sig_gen(i) = d_sig_gen(i)*dtimesub 3367 d_sig_death(i) = d_sig_death(i)*dtimesub 3368 d_sig_col(i) = d_sig_col(i)*dtimesub 3369 d_sig_spread(i) = d_sig_spread(i)*dtimesub 3370 d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) 3371 3372 3373 d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i)) 3374 !! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) 3375 !! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) 3376 d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) 3377 d_sigmaw(i) = d_sigmaw_targ 3378 !! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) 3379 3380 3381 d_dens_gen(i) = wgen(i) 3382 d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min 3383 d_dens_col(i) = - cont_fact(i)*wdens(i)**2 3384 ! 3385 d_dens_gen(i) = d_dens_gen(i)*dtimesub 3386 d_dens_death(i) = d_dens_death(i)*dtimesub 3387 d_dens_col(i) = d_dens_col(i)*dtimesub 3388 d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i) 3389 3390 3391 d_adens_death(i) = -awdens(i)/tau_prime(i) 3392 d_adens_icol(i) = cont_fact(i)*(wdens(i)-awdens(i))**2 3393 d_adens_acol(i) = - cont_fact(i)*awdens(i)**2 3394 ! 3395 d_adens_death(i) = d_adens_death(i)*dtimesub 3396 d_adens_icol(i) = d_adens_icol(i)*dtimesub 3397 d_adens_acol(i) = d_adens_acol(i)*dtimesub 3398 d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i) 3399 3400 !! 3401 d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i)) 3402 !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) 3403 d_dens_bnd(i) = d_wdens_targ - d_wdens(i) 3404 d_wdens(i) = d_wdens_targ 3405 3406 d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i)) 3407 !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) 3408 d_adens_bnd(i) = d_wdens_targ - d_awdens(i) 3409 d_awdens(i) = d_wdens_targ 3410 3411 3412 3413 ENDIF 3414 ENDDO 3415 3416 IF (prt_level >= 10) THEN 3417 print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) ', & 3418 cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) 3419 print *,'wake, wdens(1), awdens(1), d_awdens(1) ', & 3420 wdens(1), awdens(1), d_awdens(1) 3421 print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & 3422 d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) 3423 ENDIF 3424 sigmaw=sigmaw+d_sigmaw 3425 wdens=wdens+d_wdens 3426 awdens=awdens+d_awdens 3427 3428 RETURN 3429 END SUBROUTINE wake_popdyn_2 3430 3431 SUBROUTINE wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, & 3432 wdensmin, & 3433 sigmaw, asigmaw, wdens, awdens, & !! state variables 3434 gfl, agfl, cstar, cin, wape, & 3435 rad_wk, arad_wk, irad_wk, & 3436 d_sigmaw, d_asigmaw, d_wdens, d_awdens, & !! tendencies 3437 d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & 3438 d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, & 3439 d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & 3440 d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd ) 3441 3442 3443 3444 USE lmdz_wake_ini , ONLY : wake_ini 3445 USE lmdz_wake_ini , ONLY : prt_level,RG 3446 USE lmdz_wake_ini , ONLY : stark, wdens_ref 3447 USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0 3448 !! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin 3449 USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn 3450 USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max 3451 USE lmdz_wake_ini , ONLY : smallestreal 3452 3453 IMPLICIT NONE 3454 3455 INTEGER, INTENT(IN) :: klon,klev 3456 LOGICAL, INTENT(IN) :: phys_sub 3457 LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv 3458 REAL, INTENT(IN) :: dtimesub 3459 REAL, INTENT(IN) :: wdensmin 3460 REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes 3461 REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes 3462 REAL, DIMENSION (klon), INTENT(INOUT) :: asigmaw !! sigma = fractional area of active wakes 3463 REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area 3464 REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area 3465 REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes 3466 REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre 3467 3468 ! 3469 REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = mean wake radius 3470 REAL, DIMENSION (klon), INTENT(OUT) :: arad_wk !! r_A = wake radius of active wakes 3471 REAL, DIMENSION (klon), INTENT(OUT) :: irad_wk !! r_I = wake radius of inactive wakes 3472 REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front length per unit area 3473 REAL, DIMENSION (klon), INTENT(OUT) :: agfl !! LgA = gust front length of active wakes 3474 !! per unit area 3475 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_asigmaw, d_wdens, d_awdens 3476 ! Some components of the tendencies of state variables 3477 REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd 3478 REAL, DIMENSION (klon), INTENT(OUT) :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd 3479 REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd 3480 REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_acol, d_adens_icol, d_adens_bnd 3481 3482 3483 !! internal variables 3484 3485 INTEGER :: i, k 3486 REAL, DIMENSION (klon) :: iwdens, isigmaw !! inactive wake density and fractional area 3487 !! REAL, DIMENSION (klon) :: d_arad, d_irad 3488 REAL, DIMENSION (klon) :: igfl !! LgI = gust front length of inactive wakes 3489 !! per unit area 3490 REAL, DIMENSION (klon) :: s_wk !! mean area of individual wakes 3491 REAL, DIMENSION (klon) :: as_wk !! mean area of individual active wakes 3492 REAL, DIMENSION (klon) :: is_wk !! mean area of individual inactive wakes 3493 REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes 3494 REAL :: tau_wk_inv_min 3495 REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes 3496 REAL :: d_wdens_targ, d_sigmaw_targ 3497 3498 3499 !! Equations 3500 !! --------- 3501 !! Gust fronts: 3502 !! Lg_A = 2 pi r_A A 3503 !! Lg_I = 2 pi r_I I 3504 !! Lg = 2 pi r D 3505 !! 3506 !! Areas: 3507 !! s = pi r^2 3508 !! s_A = pi r_A^2 3509 !! s_I = pi r_I^2 3510 !! 3511 !! Life expectancy: 3512 !! tau_I = 3 C* ((C*/C*t)^3/2 - 1) / r_I 3513 !! 3514 !! Time deratives: 3515 !! dD/dt = B - (D-A)/tau_I - 2 Lg C* D 3516 !! dA/dt = B - A/tau_A + 2 Lg_I C* (D-A) - 2 Lg_A C* A 3517 !! dsigma/dt = B a0 - sigma_I/tau_I + Lg C* - 2 Lg_I C* (D-A) (2 s_I - a0) 3518 !! dsigma_A/dt = B a0 - sigma_A/tau_A + Lg_A C* + (Lg_A I + Lg_I A) C* s_I + 2 Lg_I C* I a0 3519 !! 3520 3521 DO i = 1, klon 3522 IF (wk_adv(i)) THEN 3523 iwdens(i) = wdens(i) - awdens(i) 3524 isigmaw(i) = sigmaw(i) - asigmaw(i) 3525 ! 3526 arad_wk(i) = max( sqrt(asigmaw(i)/(3.14*awdens(i))) , rzero) 3527 irad_wk(i) = max( sqrt((sigmaw(i)-asigmaw(i))/ & 3528 (3.14*max(smallestreal,(wdens(i)-awdens(i))))), rzero) 3529 rad_wk(i) = (awdens(i)*arad_wk(i)+(wdens(i)-awdens(i))*irad_wk(i))/wdens(i) 3530 ! 3531 s_wk(i) = 3.14*rad_wk(i)**2 3532 as_wk(i) = 3.14*arad_wk(i)**2 3533 is_wk(i) = 3.14*irad_wk(i)**2 3534 ! 3535 gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) 3536 agfl(i) = 2.*sqrt(3.14*awdens(i)*asigmaw(i)) 3537 igfl(i) = gfl(i) - agfl(i) 3538 ENDIF 3539 ENDDO 3540 3541 3542 DO i = 1, klon 3543 IF (wk_adv(i)) THEN 3544 tau_wk_inv(i) = max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.) 3545 tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub) 3546 tau_prime(i) = tau_cv 3547 3548 d_sig_gen(i) = wgen(i)*aa0 3549 d_sig_death(i) = - isigmaw(i)*tau_wk_inv_min 3550 d_sig_col(i) = - 2.*igfl(i)*cstar(i)*iwdens(i)*(2.*is_wk(i)-aa0) 3551 d_sig_spread(i) = gfl(i)*cstar(i) 3552 ! 3553 d_sig_gen(i) = d_sig_gen(i)*dtimesub 3554 d_sig_death(i) = d_sig_death(i)*dtimesub 3555 d_sig_col(i) = d_sig_col(i)*dtimesub 3556 d_sig_spread(i) = d_sig_spread(i)*dtimesub 3557 d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) 3558 #ifdef IOPHYS_WK 3559 IF (phys_sub) call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw) 3560 #endif 3561 3562 3563 d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i)) 3564 !! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) 3565 !! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) 3566 d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) 3567 d_sigmaw(i) = d_sigmaw_targ 3568 !! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) 3569 #ifdef IOPHYS_WK 3570 IF (phys_sub) THEN 3571 call iophys_ecrit('tauwk_inv',1,'tau_wk_inv_min','',tau_wk_inv_min) 3572 call iophys_ecrit('d_sigmaw',1,'d_sigmaw','',d_sigmaw) 3573 call iophys_ecrit('d_sig_gen',1,'d_sig_gen','',d_sig_gen) 3574 call iophys_ecrit('d_sig_death',1,'d_sig_death','',d_sig_death) 3575 call iophys_ecrit('d_sig_col',1,'d_sig_col','',d_sig_col) 3576 call iophys_ecrit('d_sig_spread',1,'d_sig_spread','',d_sig_spread) 3577 call iophys_ecrit('d_sig_bnd',1,'d_sig_bnd','',d_sig_bnd) 3578 ENDIF 3579 #endif 3580 d_asig_death(i) = - asigmaw(i)/tau_prime(i) 3581 d_asig_aicol(i) = (agfl(i)*iwdens(i) + igfl(i)*awdens(i))*cstar(i)*is_wk(i) 3582 d_asig_iicol(i) = 2.*igfl(i)*cstar(i)*iwdens(i)*aa0 3583 d_asig_spread(i) = agfl(i)*cstar(i) 3584 ! 3585 d_asig_death(i) = d_asig_death(i)*dtimesub 3586 d_asig_aicol(i) = d_asig_aicol(i)*dtimesub 3587 d_asig_iicol(i) = d_asig_iicol(i)*dtimesub 3588 d_asig_spread(i) = d_asig_spread(i)*dtimesub 3589 d_asigmaw(i) = d_sig_gen(i) + d_asig_death(i) + d_asig_aicol(i) + d_asig_iicol(i) + d_asig_spread(i) 3590 #ifdef IOPHYS_WK 3591 IF (phys_sub) call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw) 3592 #endif 3593 3594 d_sigmaw_targ = min(max(d_asigmaw(i),-asigmaw(i)), sigmaw(i)-asigmaw(i)) 3595 !! d_dens_bnd(i) = d_dens_bnd(i) + d_sigmaw_targ - d_sigmaw(i) 3596 d_asig_bnd(i) = d_sigmaw_targ - d_asigmaw(i) 3597 d_asigmaw(i) = d_sigmaw_targ 3598 #ifdef IOPHYS_WK 3599 IF (phys_sub) THEN 3600 call iophys_ecrit('d_asigmaw',1,'d_asigmaw','',d_asigmaw) 3601 call iophys_ecrit('d_asig_death',1,'d_asig_death','',d_asig_death) 3602 call iophys_ecrit('d_asig_aicol',1,'d_asig_aicol','',d_asig_aicol) 3603 call iophys_ecrit('d_asig_iicol',1,'d_asig_iicol','',d_asig_iicol) 3604 call iophys_ecrit('d_asig_spread',1,'d_asig_spread','',d_asig_spread) 3605 call iophys_ecrit('d_asig_bnd',1,'d_asig_bnd','',d_asig_bnd) 3606 ENDIF 3607 #endif 3608 d_dens_gen(i) = wgen(i) 3609 d_dens_death(i) = - iwdens(i)*tau_wk_inv_min 3610 d_dens_col(i) = - 2.*gfl(i)*cstar(i)*wdens(i) 3611 ! 3612 d_dens_gen(i) = d_dens_gen(i)*dtimesub 3613 d_dens_death(i) = d_dens_death(i)*dtimesub 3614 d_dens_col(i) = d_dens_col(i)*dtimesub 3615 d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i) 3616 !! 3617 d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i)) 3618 !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) 3619 d_dens_bnd(i) = d_wdens_targ - d_wdens(i) 3620 d_wdens(i) = d_wdens_targ 3621 #ifdef IOPHYS_WK 3622 IF (phys_sub) THEN 3623 call iophys_ecrit('d_wdens',1,'d_wdens','',d_wdens) 3624 call iophys_ecrit('d_dens_gen',1,'d_dens_gen','',d_dens_gen) 3625 call iophys_ecrit('d_dens_death',1,'d_dens_death','',d_dens_death) 3626 call iophys_ecrit('d_dens_col',1,'d_dens_col','',d_dens_col) 3627 ENDIF 3628 #endif 3629 3630 d_adens_death(i) = -awdens(i)/tau_prime(i) 3631 d_adens_icol(i) = 2.*igfl(i)*cstar(i)*iwdens(i) 3632 d_adens_acol(i) = - 2.*agfl(i)*cstar(i)*awdens(i) 3633 ! 3634 d_adens_death(i) = d_adens_death(i)*dtimesub 3635 d_adens_icol(i) = d_adens_icol(i)*dtimesub 3636 d_adens_acol(i) = d_adens_acol(i)*dtimesub 3637 d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i) 3638 #ifdef IOPHYS_WK 3639 IF (phys_sub) THEN 3640 call iophys_ecrit('d_awdens',1,'d_awdens','',d_awdens) 3641 call iophys_ecrit('d_adens_death',1,'d_adens_death','',d_adens_death) 3642 call iophys_ecrit('d_adens_icol',1,'d_adens_icol','',d_adens_icol) 3643 call iophys_ecrit('d_adens_acol',1,'d_adens_acol','',d_adens_acol) 3644 ENDIF 3645 #endif 3646 d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i)) 3647 !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) 3648 d_adens_bnd(i) = d_wdens_targ - d_awdens(i) 3649 d_awdens(i) = d_wdens_targ 3650 3651 !! d_irad(i) = (d_sigmaw(i)-d_asigmaw(i)-isigmaw(i)*(d_wdens(i)-awdens(i))/iwdens(i)) / & 3652 !! max(smallestreal,(2.*3.14*iwdens(i)*irad_wk(i))) 3653 !! d_arad(i) = (d_asigmaw(i)-asigmaw(i)*d_awdens(i)/awdens(i)) / & 3654 !! max(smallestreal,(2.*3.14*awdens(i)*arad_wk(i))) 3655 !! d_irad(i) = d_irad(i)*dtimesub 3656 !! d_arad(i) = d_arad(i)*dtimesub 3657 !! call iophys_ecrit('d_irad',1,'d_irad','m',d_irad) 3658 !! call iophys_ecrit('d_airad',1,'d_arad','m',d_arad) 3659 !! 3660 ENDIF 3661 ENDDO 3662 3663 IF (prt_level >= 10) THEN 3664 print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1) ', & 3665 cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1) 3666 print *,'wake, wdens(1), awdens(1), d_awdens(1) ', & 3667 wdens(1), awdens(1), d_awdens(1) 3668 print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & 3669 d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) 3670 ENDIF 3671 sigmaw=sigmaw+d_sigmaw 3672 asigmaw=asigmaw+d_asigmaw 3673 wdens=wdens+d_wdens 3674 awdens=awdens+d_awdens 3675 3676 RETURN 3677 END SUBROUTINE wake_popdyn_3 3678 2629 3679 END MODULE lmdz_wake -
LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
r5253 r5255 329 329 d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection 330 330 ! tendencies of wake fractional area and wake number per unit area: 331 d_s_wk, d_dens_a_wk, d_dens_wk, & ! due to wakes331 d_s_wk, d_s_a_wk, d_dens_a_wk, d_dens_wk, & ! due to wakes 332 332 !!! d_s_vdf, d_dens_a_vdf, d_dens_vdf, & ! due to vertical diffusion 333 333 !!! d_s_the, d_dens_a_the, d_dens_the, & ! due to thermals … … 1414 1414 endif ! 1415 1415 !======================================================================! 1416 #ifdef ISO 1417 #ifdef ISOVERIF 1418 write(*,*) 'physiq 1421b: qx(1,1,:)=',qx(1,1,:) 1419 #endif 1420 #endif 1416 1421 1417 1422 … … 4445 4450 dt_a, dq_a, cv_gen, & 4446 4451 sigd, cin, & 4447 wake_deltat, wake_deltaq, wake_s, awake_ dens,wake_dens, &4452 wake_deltat, wake_deltaq, wake_s, awake_s, wake_dens, awake_dens, & 4448 4453 wake_dth, wake_h, & 4449 4454 !! wake_pe, wake_fip, wake_gfl, & … … 4455 4460 wake_omg, wake_dp_deltomg, & 4456 4461 wake_spread, wake_Cstar, d_deltat_wk_gw, & 4457 d_deltat_wk, d_deltaq_wk, d_s_wk, d_ dens_a_wk, d_dens_wk &4462 d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk & 4458 4463 #ifdef ISO 4459 4464 & ,xt_seri,dxt_dwn,dxt_a & … … 4495 4500 4496 4501 CALL add_wake_tend & 4497 (d_deltat_wk, d_deltaq_wk, d_s_wk, d_ dens_a_wk, d_dens_wk, wake_k, &4502 (d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk, wake_k, & 4498 4503 'wake', abortphy & 4499 4504 #ifdef ISO … … 4770 4775 ,zqla,ztva & 4771 4776 #ifdef ISO 4772 & ,xt_ env,d_xt_ajs &4777 & ,xt_therm,d_xt_ajs & 4773 4778 #ifdef DIAGISO 4774 4779 & ,q_the,xt_the & … … 6557 6562 IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',& 6558 6563 minval(d_q_emiss),maxval(d_q_emiss) 6564 #ifdef ISO 6565 abort_message='isotopes pas encore dans emission volc H2O' 6566 CALL abort_physic(modname,abort_message,1) 6567 #endif 6559 6568 6560 6569 CALL add_phys_tend(du0, dv0, dt0, d_q_emiss, dql0, dqi0, dqbs0, paprs, & … … 6948 6957 6949 6958 print*,'avt add phys rep',abortphy 6950 6959 #ifdef ISO 6960 abort_message='isotopes pas encore dans rep' 6961 CALL abort_physic(modname,abort_message,1) 6962 #endif 6951 6963 CALL add_phys_tend & 6952 6964 (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,dqbs0,paprs,&
Note: See TracChangeset
for help on using the changeset viewer.
