Changeset 5256 for LMDZ6/trunk/libf/phylmdiso
- Timestamp:
- Oct 22, 2024, 4:06:28 PM (5 weeks ago)
- Location:
- LMDZ6/trunk/libf/phylmdiso
- Files:
-
- 4 edited
-
calwake.F90 (modified) (13 diffs)
-
isotopes_routines_mod.F90 (modified) (2 diffs, 1 prop)
-
lmdz_wake.F90 (modified) (96 diffs)
-
physiq_mod.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmdiso/calwake.F90
r5255 r5256 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_ s, wake_dens, awake_dens, &8 wake_deltat, wake_deltaq, wake_s, awake_dens, wake_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 s, wake_ddens, awake_ddens &16 wake_ddeltat, wake_ddeltaq, wake_ds, awake_ddens, wake_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 , awake_s70 REAL, DIMENSION(klon), INTENT (INOUT) :: wake_dens, awake_dens69 REAL, DIMENSION(klon), INTENT (INOUT) :: wake_s 70 REAL, DIMENSION(klon), INTENT (INOUT) :: awake_dens, wake_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 s, wake_ddens, awake_ddens90 REAL, DIMENSION(klon), INTENT (OUT) :: wake_ds, awake_ddens, wake_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 sigmaw, wdens, awdens117 REAL, DIMENSION(klon) :: sigmaw, awdens, wdens 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 sigmaw, d_wdens, d_awdens128 REAL, DIMENSION(klon) :: d_sigmaw, d_awdens, d_wdens 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, awake_s, wgen input ', wake_s(1), awake_s(1), wgen(1)144 print *, '-> calwake, wake_s, wgen input ', wake_s(1), wgen(1) 145 145 ENDIF 146 146 … … 191 191 d_deltaqw(:,:) = 0. 192 192 d_sigmaw(:) = 0. 193 d_a sigmaw(:) = 0.193 d_awdens(:) = 0. 194 194 d_wdens(:) = 0. 195 d_awdens(:) = 0.196 195 #ifdef ISO 197 196 dxtls(:,:,:) = 0. … … 228 227 229 228 DO i = 1, klon 230 sigmaw(i) 231 asigmaw(i) = awake_s(i)232 END DO 233 234 DO i = 1, klon229 sigmaw(i) = wake_s(i) 230 END DO 231 232 DO i = 1, klon 233 awdens(i) = max(0., awake_dens(i)) 235 234 wdens(i) = max(0., wake_dens(i)) 236 awdens(i) = max(0., awake_dens(i))237 235 END DO 238 236 … … 274 272 #endif 275 273 276 277 CALL wake(klon,klev,znatsurf, p, ph, pi, dtime, & 274 CALL wake(znatsurf, p, ph, pi, dtime, & 278 275 te, qe, omgbe, & 279 276 dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, & 280 277 sigd0, Cin, & 281 dtw, dqw, sigmaw, a sigmaw, wdens, awdens, &! state variables278 dtw, dqw, sigmaw, awdens, wdens, & ! state variables 282 279 dth, hw, wape, fip, gfl, & 283 280 dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, & 284 281 dtke, dqke, omg, dp_deltomg, spread, cstar, & 285 282 d_deltat_gw, & 286 d_deltatw, d_deltaqw, d_sigmaw, d_a sigmaw, d_wdens, d_awdens &! tendencies283 d_deltatw, d_deltaqw, d_sigmaw, d_awdens, d_wdens & ! tendencies 287 284 #ifdef ISO 288 285 , xte,dxtdwn,dxta,dxtw & … … 389 386 IF (ktopw(i)>0) THEN 390 387 wake_ds(i) = d_sigmaw(i)*dtime 391 awake_ds(i) = d_asigmaw(i)*dtime392 388 awake_ddens(i) = d_awdens(i)*dtime 393 389 wake_ddens(i) = d_wdens(i)*dtime 394 390 ELSE 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) 391 wake_ds(i) = -wake_s(i) 392 wake_ddens(i)= -wake_dens(i) 399 393 END IF 400 394 END DO … … 427 421 DO i = 1, klon 428 422 wake_s(i) = sigmaw(i) 429 awake_s(i) = asigmaw(i)430 423 awake_dens(i) = awdens(i) 431 424 wake_dens(i) = wdens(i) … … 442 435 ENDIF ! (first) 443 436 !>jyg 444 IF (prt_level >= 10) THEN445 print *, 'calwake ->, wake_s, awake_s ', wake_s(1), awake_s(1)446 ENDIF447 437 448 438 RETURN 449 439 END SUBROUTINE calwake 450 440 451 -
LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90
-
Property
svn:keywords
set to
Id
r5255 r5256 1 1 #ifdef ISO 2 ! $Id: $ 3 2 ! $Id$ 4 3 MODULE isotopes_routines_mod 5 4 USE infotrac_phy, ONLY: niso, ntraciso=>ntiso, index_trac=>itZonIso, ntraceurs_zone=>nzone … … 18822 18821 if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then 18823 18822 if (q.gt.ridicule) then 18824 !write(*,*) 'xt,q=',xt,q18825 !write(*,*) 'alpha=',alpha18826 !write(*,*) 'toce,kcin,h0=',toce,kcin,h018827 !write(*,*) 'RMerlivat=',RMerlivat18823 write(*,*) 'xt,q=',xt,q 18824 write(*,*) 'alpha=',alpha 18825 write(*,*) 'toce,kcin,h0=',toce,kcin,h0 18826 write(*,*) 'RMerlivat=',RMerlivat 18828 18827 call iso_verif_aberrant_encadre( xt/q, 'isotopes_routines_mod 18930b: iso_init_ideal') 18829 18828 endif -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmdiso/lmdz_wake.F90
r5255 r5256 4 4 5 5 CONTAINS 6 7 SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, & 8 tb0, qb0, omgb, & 6 SUBROUTINE wake(znatsurf, p, ph, pi, dtime, & 7 te0, qe0, omgb, & 9 8 dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, & 10 9 sigd_con, Cin, & 11 deltatw, deltaqw, sigmaw, a sigmaw, wdens, awdens, & ! state variables10 deltatw, deltaqw, sigmaw, awdens, wdens, & ! state variables 12 11 dth, hw, wape, fip, gfl, & 13 dtls, dqls, ktopw, omgbdth, dp_omgb, t x, qx, &14 dtke, dqke, omg, dp_deltomg, wkspread, cstar, &15 d_deltat_gw, & ! tendencies16 d_deltatw2, d_deltaqw2, d_sigmaw2, d_a sigmaw2, d_wdens2, d_awdens2 & ! tendencies17 18 #ifdef ISO 19 ,xt b0,dxtdwn,dxta,deltaxtw &20 ,dxtls,xt x,dxtke,d_deltaxtw2 &12 dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, & 13 dtke, dqke, omg, dp_deltomg, spread, cstar, & 14 d_deltat_gw, & 15 d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2 & ! tendencies 16 17 #ifdef ISO 18 ,xte0,dxtdwn,dxta,deltaxtw & 19 ,dxtls,xtu,dxtke,d_deltaxtw2 & 21 20 #endif 22 21 ) … … 32 31 ! ************************************************************** 33 32 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 33 USE ioipsl_getin_p_mod, ONLY : getin_p 34 USE dimphy 35 use mod_phys_lmdz_para 36 USE print_control_mod, ONLY: prt_level 46 37 #ifdef ISO 47 38 USE infotrac_phy, ONLY : ntraciso=>ntiso … … 64 55 ! deltaqw : specific humidity difference between wake and off-wake regions 65 56 ! sigmaw : fractional area covered by wakes. 66 ! asigmaw : fractional area covered by active wakes.67 57 ! wdens : number of wakes per unit area 68 ! awdens : number of active wakes per unit area69 58 70 59 ! Variable de sortie : … … 84 73 ! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) 85 74 ! dp_deltomg : vertical gradient of omg (s-1) 86 ! wkspread : spreading term in d_t_wake and d_q_wake75 ! spread : spreading term in d_t_wake and d_q_wake 87 76 ! deltatw : updated temperature difference (T_w-T_u). 88 77 ! deltaqw : updated humidity difference (q_w-q_u). 89 78 ! sigmaw : updated wake fractional area. 90 ! asigmaw : updated active wake fractional area.91 79 ! d_deltat_gw : delta T tendency due to GW 92 80 … … 94 82 95 83 ! aire : aire de la maille 96 ! t b0 : horizontal average of temperature(K)97 ! q b0 : horizontal average of humidity(kg/kg)84 ! te0 : temperature dans l'environnement (K) 85 ! qe0 : humidite dans l'environnement (kg/kg) 98 86 ! omgb : vitesse verticale moyenne sur la maille (Pa/s) 99 87 ! dtdwn: source de chaleur due aux descentes (K/s) … … 115 103 ! Variables internes : 116 104 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 105 ! rhow : masse volumique de la poche froide 106 ! rho : environment density at P levels 107 ! rhoh : environment density at Ph levels 108 ! te : environment temperature | may change within 109 ! qe : environment humidity | sub-time-stepping 110 ! the : environment potential temperature 111 ! thu : potential temperature in undisturbed area 112 ! tu : temperature in undisturbed area 113 ! qu : humidity in undisturbed area 125 114 ! dp_omgb: vertical gradient og LS omega 126 115 ! omgbw : wake average vertical omega … … 128 117 ! omgbdq : flux of Delta_q transported by LS omega 129 118 ! dth : potential temperature diff. wake-undist. 130 ! th1 : first pot. temp. for vertical advection (=th x)119 ! th1 : first pot. temp. for vertical advection (=thu) 131 120 ! th2 : second pot. temp. for vertical advection (=thw) 132 121 ! q1 : first humidity for vertical advection 133 122 ! q2 : second humidity for vertical advection 134 ! d_deltatw : redistribution term for deltatw135 ! d_deltaqw : redistribution term for deltaqw136 ! deltatw0 : initial deltatw137 ! deltaqw0 : initial deltaqw123 ! d_deltatw : terme de redistribution pour deltatw 124 ! d_deltaqw : terme de redistribution pour deltaqw 125 ! deltatw0 : deltatw initial 126 ! deltaqw0 : deltaqw initial 138 127 ! hw0 : wake top hight (defined as the altitude at which deltatw=0) 139 128 ! amflux : horizontal mass flux through wake boundary 140 129 ! wdens_ref: initial number of wakes per unit area (3D) or per 141 ! 142 ! Tgw : 1 sur la p eriode de onde de gravite143 ! Cgw : vitesse de propagation de onde de gravit e144 ! LL : distance between 2 wakes130 ! unit length (2D), at the beginning of each time step 131 ! Tgw : 1 sur la période de onde de gravité 132 ! Cgw : vitesse de propagation de onde de gravité 133 ! LL : distance entre 2 poches 145 134 146 135 ! ------------------------------------------------------------------------- 147 ! D eclaration de variables136 ! Déclaration de variables 148 137 ! ------------------------------------------------------------------------- 149 138 139 include "YOMCST.h" 140 include "cvthermo.h" 150 141 151 142 ! Arguments en entree 152 143 ! -------------------- 153 144 154 INTEGER, INTENT(IN) :: klon,klev155 145 INTEGER, DIMENSION (klon), INTENT(IN) :: znatsurf 156 146 REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi … … 158 148 REAL, DIMENSION (klon, klev), INTENT(IN) :: omgb 159 149 REAL, INTENT(IN) :: dtime 160 REAL, DIMENSION (klon, klev), INTENT(IN) :: t b0, qb0150 REAL, DIMENSION (klon, klev), INTENT(IN) :: te0, qe0 161 151 REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn 162 152 REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup … … 166 156 REAL, DIMENSION (klon), INTENT(IN) :: Cin 167 157 #ifdef ISO 168 REAL, DIMENSION (ntraciso,klon, klev), INTENT(IN) :: xt b0158 REAL, DIMENSION (ntraciso,klon, klev), INTENT(IN) :: xte0 169 159 REAL, DIMENSION (ntraciso,klon, klev), INTENT(IN) :: dxtdwn 170 160 REAL, DIMENSION (ntraciso,klon, klev), INTENT(IN) :: dxta … … 176 166 REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw 177 167 REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw 178 REAL, DIMENSION (klon), INTENT(INOUT) :: a sigmaw168 REAL, DIMENSION (klon), INTENT(INOUT) :: awdens 179 169 REAL, DIMENSION (klon), INTENT(INOUT) :: wdens 180 REAL, DIMENSION (klon), INTENT(INOUT) :: awdens181 170 #ifdef ISO 182 171 REAL, DIMENSION (ntraciso,klon, klev), INTENT(INOUT) :: deltaxtw … … 187 176 188 177 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth 189 REAL, DIMENSION (klon, klev), INTENT(OUT) :: t x, qx178 REAL, DIMENSION (klon, klev), INTENT(OUT) :: tu, qu 190 179 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls 191 180 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke 192 REAL, DIMENSION (klon, klev), INTENT(OUT) :: wkspread ! unused (jyg)181 REAL, DIMENSION (klon, klev), INTENT(OUT) :: spread ! unused (jyg) 193 182 REAL, DIMENSION (klon, klev), INTENT(OUT) :: omgbdth, omg 194 183 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg 184 REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw 195 185 REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar 196 186 INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw 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 187 ! Tendencies of state variables 200 188 REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2 201 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2 202 203 #ifdef ISO 204 REAL, DIMENSION (ntraciso,klon, klev), INTENT(OUT) :: xtx 189 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_awdens2, d_wdens2 190 #ifdef ISO 191 REAL, DIMENSION (ntraciso,klon, klev), INTENT(OUT) :: xtu 205 192 REAL, DIMENSION (ntraciso,klon, klev), INTENT(OUT) :: dxtls 206 193 REAL, DIMENSION (ntraciso,klon, klev), INTENT(OUT) :: dxtke … … 211 198 ! ------------------- 212 199 213 ! Variables a fixer 200 ! Variables à fixer 201 INTEGER, SAVE :: igout 202 !$OMP THREADPRIVATE(igout) 203 LOGICAL, SAVE :: first = .TRUE. 204 !$OMP THREADPRIVATE(first) 205 !jyg< 206 !! REAL, SAVE :: stark, wdens_ref, coefgw, alpk 207 REAL, SAVE, DIMENSION(2) :: wdens_ref 208 REAL, SAVE :: stark, coefgw, alpk 209 !>jyg 210 REAL, SAVE :: crep_upper, crep_sol 211 !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol) 212 213 REAL, SAVE :: tau_cv 214 !$OMP THREADPRIVATE(tau_cv) 215 REAL, SAVE :: rzero, aa0 ! minimal wake radius and area 216 !$OMP THREADPRIVATE(rzero, aa0) 217 218 LOGICAL, SAVE :: flag_wk_check_trgl 219 !$OMP THREADPRIVATE(flag_wk_check_trgl) 220 INTEGER, SAVE :: iflag_wk_check_trgl 221 !$OMP THREADPRIVATE(iflag_wk_check_trgl) 222 INTEGER, SAVE :: iflag_wk_pop_dyn 223 !$OMP THREADPRIVATE(iflag_wk_pop_dyn) 214 224 215 225 REAL :: delta_t_min 226 INTEGER :: nsub 216 227 REAL :: dtimesub 228 REAL, SAVE :: wdensmin 229 !$OMP THREADPRIVATE(wdensmin) 230 REAL, SAVE :: sigmad, hwmin, wapecut, cstart 231 !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart) 232 REAL, SAVE :: sigmaw_max 233 !$OMP THREADPRIVATE(sigmaw_max) 234 REAL, SAVE :: dens_rate 235 !$OMP THREADPRIVATE(dens_rate) 217 236 REAL :: wdens0 218 237 ! IM 080208 … … 222 241 REAL, DIMENSION (klon, klev) :: deltatw0 223 242 REAL, DIMENSION (klon, klev) :: deltaqw0 224 REAL, DIMENSION (klon, klev) :: tb, qb 243 REAL, DIMENSION (klon, klev) :: te, qe 244 !! REAL, DIMENSION (klon) :: sigmaw1 225 245 #ifdef ISO 226 246 REAL, DIMENSION (ntraciso,klon, klev) :: deltaxtw0 227 REAL, DIMENSION (ntraciso,klon, klev) :: xt b228 #endif 229 230 ! Variables liees a la dynamique de population 1247 REAL, DIMENSION (ntraciso,klon, klev) :: xte 248 #endif 249 250 ! Variables liees a la dynamique de population 231 251 REAL, DIMENSION(klon) :: act 232 252 REAL, DIMENSION(klon) :: rad_wk, tau_wk_inv 233 253 REAL, DIMENSION(klon) :: f_shear 234 254 REAL, DIMENSION(klon) :: drdt 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 255 REAL, DIMENSION(klon) :: d_sig_gen, d_sig_death, d_sig_col 243 256 REAL, DIMENSION(klon) :: wape1_act, wape2_act 244 257 LOGICAL, DIMENSION (klon) :: kill_wake 258 INTEGER, SAVE :: iflag_wk_act 259 !$OMP THREADPRIVATE(iflag_wk_act) 245 260 REAL :: drdt_pos 246 261 REAL :: tau_wk_inv_min 247 ! Some components of the tendencies of state variables248 REAL, DIMENSION (klon) :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2249 REAL, DIMENSION (klon) :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2250 REAL, DIMENSION (klon) :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2251 REAL, DIMENSION (klon) :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2252 262 253 263 ! Variables pour les GW … … 258 268 259 269 ! Variables liees au calcul de hw 260 REAL, DIMENSION (klon) :: ptop 270 REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new 261 271 REAL, DIMENSION (klon) :: sum_dth 262 272 REAL, DIMENSION (klon) :: dthmin … … 270 280 ! Sub-timestep tendencies and related variables 271 281 REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw 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 282 REAL, DIMENSION (klon, klev) :: d_te, d_qe 283 REAL, DIMENSION (klon) :: d_awdens, d_wdens 284 REAL, DIMENSION (klon) :: d_sigmaw, alpha 281 285 REAL, DIMENSION (klon) :: q0_min, q1_min 282 286 LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw 283 287 #ifdef ISO 284 288 REAL, DIMENSION (ntraciso,klon, klev) :: d_deltaxtw 285 REAL, DIMENSION (ntraciso,klon, klev) :: d_xtb 286 #endif 289 REAL, DIMENSION (ntraciso,klon, klev) :: d_xte 290 #endif 291 REAL, SAVE :: epsilon 292 !$OMP THREADPRIVATE(epsilon) 293 DATA epsilon/1.E-15/ 287 294 288 295 ! Autres variables internes 289 INTEGER ::isubstep, k, i, igout 290 291 REAL :: wdensmin 292 296 INTEGER ::isubstep, k, i 297 298 REAL :: wdens_targ 293 299 REAL :: sigmaw_targ 294 REAL :: wdens_targ 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 300 301 REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu 302 REAL, DIMENSION (klon) :: sum_dq, sum_rho 300 303 REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn 301 REAL, DIMENSION (klon) :: av_th x, av_tx, av_qx, av_thvx302 REAL, DIMENSION (klon) :: av_dth, av_dq 304 REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu 305 REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho 303 306 REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn 304 305 REAL, DIMENSION (klon, klev) :: rho 307 ! pas besoin d'isos ici 308 309 REAL, DIMENSION (klon, klev) :: rho, rhow 306 310 REAL, DIMENSION (klon, klev+1) :: rhoh 311 REAL, DIMENSION (klon, klev) :: rhow_moyen 307 312 REAL, DIMENSION (klon, klev) :: zh 308 313 REAL, DIMENSION (klon, klev+1) :: zhh 309 310 REAL, DIMENSION (klon, klev) :: thb, thx 314 REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2 315 316 REAL, DIMENSION (klon, klev) :: the, thu 311 317 312 318 REAL, DIMENSION (klon, klev) :: omgbw … … 343 349 REAL, DIMENSION (klon, klev) :: detr 344 350 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 REAL, DIMENSION(klon) :: sigmaw_in ! pour les prints 352 REAL, DIMENSION(klon) :: awdens_in, wdens_in ! pour les prints 364 353 365 354 ! ------------------------------------------------------------------------- 366 355 ! Initialisations 367 356 ! ------------------------------------------------------------------------- 357 358 ! print*, 'wake initialisations' 359 !#ifdef ISOVERIF 360 ! write(*,*) 'wake 358: entree' 361 !#endif 362 363 ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10. 364 ! ------------------------------------------------------------------------- 365 366 !! DATA wapecut, sigmad, hwmin/5., .02, 10./ 367 !! DATA wapecut, sigmad, hwmin/1., .02, 10./ 368 DATA sigmad, hwmin/.02, 10./ 369 !! DATA wdensmin/1.e-12/ 370 DATA wdensmin/1.e-14/ 371 ! cc nrlmd 372 DATA sigmaw_max/0.4/ 373 DATA dens_rate/0.1/ 374 ! cc 375 ! Longueur de maille (en m) 376 ! ------------------------------------------------------------------------- 377 368 378 ! ALON = 3.e5 369 379 ! alon = 1.E6 … … 372 382 f_shear(:) = 1. ! 0. for strong shear, 1. for weak shear 373 383 384 374 385 ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) 375 386 376 ! coefgw : Coefficient pour les ondes de gravit e387 ! coefgw : Coefficient pour les ondes de gravité 377 388 ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) 378 ! wdens : Densit esurfacique de poche froide389 ! wdens : Densité surfacique de poche froide 379 390 ! ------------------------------------------------------------------------- 380 391 … … 389 400 ! alpk = 0.5 390 401 ! alpk = 0.05 402 403 if (first) then 404 405 igout = klon/2+1/klon 406 407 crep_upper = 0.9 408 crep_sol = 1.0 409 410 ! Get wapecut from parameter file 411 wapecut = 1. 412 CALL getin_p('wapecut', wapecut) 413 414 ! cc nrlmd Lecture du fichier wake_param.data 415 stark=0.33 416 CALL getin_p('stark',stark) 417 cstart = stark*sqrt(2.*wapecut) 418 419 alpk=0.25 420 CALL getin_p('alpk',alpk) 421 !jyg< 422 !! wdens_ref=8.E-12 423 !! CALL getin_p('wdens_ref',wdens_ref) 424 wdens_ref(1)=8.E-12 425 wdens_ref(2)=8.E-12 426 CALL getin_p('wdens_ref_o',wdens_ref(1)) !wake number per unit area ; ocean 427 CALL getin_p('wdens_ref_l',wdens_ref(2)) !wake number per unit area ; land 428 !>jyg 391 429 ! 392 igout = klon/2+1/klon 430 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 431 !!!!!!!!! Population dynamics parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!! 432 !------------------------------------------------------------------------ 433 434 iflag_wk_pop_dyn = 0 435 CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed 436 ! and wdens prognostic 437 iflag_wk_act = 0 438 CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0. 439 ! 1: act(:)=1. 440 ! 2: act(:)=f(Wape) 441 442 rzero = 5000. 443 CALL getin_p('rzero_wk', rzero) 444 aa0 = 3.14*rzero*rzero 393 445 ! 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) 446 tau_cv = 4000. 447 CALL getin_p('tau_cv', tau_cv) 448 449 !------------------------------------------------------------------------ 450 451 coefgw=4. 452 CALL getin_p('coefgw',coefgw) 453 454 WRITE(*,*) 'stark=', stark 455 WRITE(*,*) 'alpk=', alpk 456 !jyg< 457 !! WRITE(*,*) 'wdens_ref=', wdens_ref 458 WRITE(*,*) 'wdens_ref_o=', wdens_ref(1) 459 WRITE(*,*) 'wdens_ref_l=', wdens_ref(2) 460 !>jyg 461 WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn 462 WRITE(*,*) 'iflag_wk_act',iflag_wk_act 463 WRITE(*,*) 'coefgw=', coefgw 464 465 flag_wk_check_trgl=.false. 466 CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl) 467 WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl 468 WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot' 469 iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1 470 CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl) 471 WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl 472 473 first=.false. 474 endif 409 475 410 476 IF (iflag_wk_pop_dyn == 0) THEN … … 419 485 !>jyg 420 486 ENDIF ! (iflag_wk_pop_dyn == 0) 421 !422 IF (iflag_wk_pop_dyn >=1) THEN423 IF (iflag_wk_pop_dyn == 3) THEN424 wdensmin = wdensthreshold425 ELSE426 wdensmin = wdensinit427 ENDIF428 ENDIF429 487 430 488 ! print*,'stark',stark … … 436 494 ! ------------------------------------------------------------------------- 437 495 438 496 delta_t_min = 0.2 439 497 440 498 ! 1. - Save initial values, initialize tendencies, initialize output fields … … 447 505 !! deltatw0(i, k) = deltatw(i, k) 448 506 !! deltaqw0(i, k) = deltaqw(i, k) 449 !! t b(i, k) = tb0(i, k)450 !! q b(i, k) = qb0(i, k)507 !! te(i, k) = te0(i, k) 508 !! qe(i, k) = qe0(i, k) 451 509 !! dtls(i, k) = 0. 452 510 !! dqls(i, k) = 0. 453 511 !! d_deltat_gw(i, k) = 0. 454 !! d_t b(i, k) = 0.455 !! d_q b(i, k) = 0.512 !! d_te(i, k) = 0. 513 !! d_qe(i, k) = 0. 456 514 !! d_deltatw(i, k) = 0. 457 515 !! d_deltaqw(i, k) = 0. … … 465 523 deltatw0(:,:) = deltatw(:,:) 466 524 deltaqw0(:,:) = deltaqw(:,:) 467 t b(:,:) = tb0(:,:)468 q b(:,:) = qb0(:,:)525 te(:,:) = te0(:,:) 526 qe(:,:) = qe0(:,:) 469 527 dtls(:,:) = 0. 470 528 dqls(:,:) = 0. 471 529 d_deltat_gw(:,:) = 0. 472 d_t b(:,:) = 0.473 d_q b(:,:) = 0.530 d_te(:,:) = 0. 531 d_qe(:,:) = 0. 474 532 d_deltatw(:,:) = 0. 475 533 d_deltaqw(:,:) = 0. 476 534 d_deltatw2(:,:) = 0. 477 d_deltaqw2(:,:) = 0. 535 d_deltaqw2(:,:) = 0. 478 536 #ifdef ISO 479 537 deltaxtw0(:,:,:)= deltaxtw(:,:,:) 480 xt b(:,:,:) = xtb0(:,:,:)538 xte(:,:,:) = xte0(:,:,:) 481 539 dxtls(:,:,:) = 0. 482 d_xt b(:,:,:) = 0.540 d_xte(:,:,:) = 0. 483 541 d_deltaxtw(:,:,:) = 0. 484 542 d_deltaxtw2(:,:,:)=0. 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. 543 xt1(:,:,:) = 0. 544 xt2(:,:,:)=0. 545 ! init non indispensable mais facilite verif iso 546 q1(:,:)=0.0 547 q2(:,:)=0.0 548 #endif 508 549 509 550 IF (iflag_wk_act == 0) THEN … … 516 557 !! sigmaw_in(i) = sigmaw(i) 517 558 !! END DO 518 sigmaw_in(:) = sigmaw(:) 519 asigmaw_in(:) = asigmaw(:) 559 sigmaw_in(:) = sigmaw(:) 520 560 !>jyg 561 562 ! sigmaw1=sigmaw 563 ! IF (sigd_con.GT.sigmaw1) THEN 564 ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con 565 ! ENDIF 566 IF (iflag_wk_pop_dyn >=1) THEN 567 DO i = 1, klon 568 wdens_targ = max(wdens(i),wdensmin) 569 d_wdens2(i) = wdens_targ - wdens(i) 570 wdens(i) = wdens_targ 571 END DO 572 ELSE 573 DO i = 1, klon 574 d_awdens2(i) = 0. 575 d_wdens2(i) = 0. 576 END DO 577 ENDIF ! (iflag_wk_pop_dyn >=1) 578 ! 579 DO i = 1, klon 580 ! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) 581 !jyg< 582 !! sigmaw(i) = amax1(sigmaw(i), sigmad) 583 !! sigmaw(i) = amin1(sigmaw(i), 0.99) 584 sigmaw_targ = min(max(sigmaw(i), sigmad),0.99) 585 d_sigmaw2(i) = sigmaw_targ - sigmaw(i) 586 sigmaw(i) = sigmaw_targ 587 !>jyg 588 END DO 589 521 590 ! 522 591 IF (iflag_wk_pop_dyn >= 1) THEN … … 528 597 ENDIF ! (iflag_wk_pop_dyn >= 1) 529 598 530 531 ! sigmaw1=sigmaw532 ! IF (sigd_con.GT.sigmaw1) THEN533 ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con534 ! ENDIF535 IF (iflag_wk_pop_dyn >= 1) THEN536 DO i = 1, klon537 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) THEN542 !! 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_targ547 ELSE548 d_dens_bnd2(i) = 0.549 d_wdens2(i) = 0.550 ENDIF !! (wdens(i) < wdensthreshold)551 END DO552 IF (iflag_wk_pop_dyn >= 2) THEN553 DO i = 1, klon554 IF (awdens(i) < wdensthreshold) THEN555 !! 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_targ560 ELSE561 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_targ565 ENDIF566 END DO567 ENDIF ! (iflag_wk_pop_dyn >= 2)568 ELSE569 DO i = 1, klon570 d_awdens2(i) = 0.571 d_wdens2(i) = 0.572 END DO573 ENDIF ! (iflag_wk_pop_dyn >= 1)574 !575 DO i = 1, klon576 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_targ580 END DO581 !582 IF (iflag_wk_pop_dyn == 3) THEN583 DO i = 1, klon584 IF ((wdens(i)-awdens(i)) <= smallestreal) THEN585 sigmaw_targ = sigmaw(i)586 ELSE587 sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))588 ENDIF589 d_asig_bnd2(i) = sigmaw_targ - asigmaw(i)590 d_asigmaw2(i) = sigmaw_targ - asigmaw(i)591 asigmaw(i) = sigmaw_targ592 END DO593 ENDIF ! (iflag_wk_pop_dyn == 3)594 595 599 wape(:) = 0. 596 600 wape2(:) = 0. 597 601 d_sigmaw(:) = 0. 598 d_asigmaw(:) = 0.599 602 ktopw(:) = 0 600 603 ! 601 604 !<jyg 602 605 dth(:,:) = 0. 603 t x(:,:) = 0.604 q x(:,:) = 0.606 tu(:,:) = 0. 607 qu(:,:) = 0. 605 608 dtke(:,:) = 0. 606 609 dqke(:,:) = 0. 607 wkspread(:,:) = 0.610 spread(:,:) = 0. 608 611 omgbdth(:,:) = 0. 609 612 omg(:,:) = 0. … … 623 626 omgbdq(:,:) = 0. 624 627 #ifdef ISO 625 xt x(:,:,:) = 0.628 xtu(:,:,:) = 0. 626 629 dxtke(:,:,:) = 0. 627 630 omgbdxt(:,:,:) = 0. … … 652 655 ktop(i) = 0 653 656 kupper(i) = 0 654 sum_th x(i) = 0.655 sum_t x(i) = 0.656 sum_q x(i) = 0.657 sum_thv x(i) = 0.657 sum_thu(i) = 0. 658 sum_tu(i) = 0. 659 sum_qu(i) = 0. 660 sum_thvu(i) = 0. 658 661 sum_dth(i) = 0. 659 662 sum_dq(i) = 0. 663 sum_rho(i) = 0. 660 664 sum_dtdwn(i) = 0. 661 665 sum_dqdwn(i) = 0. 662 666 663 av_th x(i) = 0.664 av_t x(i) = 0.665 av_q x(i) = 0.666 av_thv x(i) = 0.667 av_thu(i) = 0. 668 av_tu(i) = 0. 669 av_qu(i) = 0. 670 av_thvu(i) = 0. 667 671 av_dth(i) = 0. 668 672 av_dq(i) = 0. 673 av_rho(i) = 0. 669 674 av_dtdwn(i) = 0. 670 675 av_dqdwn(i) = 0. … … 672 677 END DO 673 678 679 680 #ifdef ISOVERIF 681 ! TODO 682 #endif 683 674 684 ! Distance between wakes 675 685 DO i = 1, klon … … 680 690 DO k = 1, klev 681 691 DO i = 1, klon 682 ! write(*,*)'wake 1',i,k, RD,tb(i,k)683 rho(i, k) = p(i, k)/( RD*tb(i,k))692 ! write(*,*)'wake 1',i,k,rd,te(i,k) 693 rho(i, k) = p(i, k)/(rd*te(i,k)) 684 694 ! write(*,*)'wake 2',rho(i,k) 685 695 IF (k==1) THEN 686 ! write(*,*)'wake 3',i,k,rd,t b(i,k)687 rhoh(i, k) = ph(i, k)/( RD*tb(i,k))688 ! write(*,*)'wake 4',i,k,rd,t b(i,k)696 ! write(*,*)'wake 3',i,k,rd,te(i,k) 697 rhoh(i, k) = ph(i, k)/(rd*te(i,k)) 698 ! write(*,*)'wake 4',i,k,rd,te(i,k) 689 699 zhh(i, k) = 0 690 700 ELSE 691 ! write(*,*)'wake 5',rd,(t b(i,k)+tb(i,k-1))692 rhoh(i, k) = ph(i, k)*2./( RD*(tb(i,k)+tb(i,k-1)))701 ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1)) 702 rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) 693 703 ! write(*,*)'wake 6',(-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)704 zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) 695 705 END IF 696 706 ! write(*,*)'wake 7',ppi(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))) 707 the(i, k) = te(i, k)/ppi(i, k) 708 thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 709 tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) 710 qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) 711 ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k))) 712 rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) 702 713 dth(i, k) = deltatw(i, k)/ppi(i, k) 703 714 #ifdef ISO 704 715 do ixt=1,ntraciso 705 xt x(ixt,i,k) = xtb(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)716 xtu(ixt,i,k) = xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i) 706 717 enddo !do ixt=1,ntraciso 707 718 #ifdef ISOVERIF 708 719 if (iso_eau.gt.0) then 709 720 call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 723a') 710 call iso_verif_egalite(q x(i,k),xtx(iso_eau,i,k),'wake 723b')721 call iso_verif_egalite(qu(i,k),xtu(iso_eau,i,k),'wake 723b') 711 722 endif 712 723 if (iso_HDO.gt.0) then 713 if (iso_verif_aberrant_enc_choix_nostop(xt x(iso_hdo,i,k),qx(i,k),ridicule,deltalim, &714 'wake 723c xt x').eq.1) then724 if (iso_verif_aberrant_enc_choix_nostop(xtu(iso_hdo,i,k),qu(i,k),ridicule,deltalim, & 725 'wake 723c xtu').eq.1) then 715 726 stop 716 727 endif … … 721 732 END DO 722 733 734 735 723 736 DO k = 1, klev - 1 724 737 DO i = 1, klon … … 726 739 n2(i, k) = 0 727 740 ELSE 728 n2(i, k) = amax1(0., - RG**2/thb(i,k)*rho(i,k)*(thb(i,k+1)-thb(i,k-1))/ &741 n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ & 729 742 (p(i,k+1)-p(i,k-1))) 730 743 END IF … … 743 756 END DO 744 757 745 758 ! Calcul de la masse volumique moyenne de la colonne (bdlmd) 759 ! ----------------------------------------------------------------- 760 761 DO k = 1, klev 762 DO i = 1, klon 763 epaisseur1(i, k) = 0. 764 epaisseur2(i, k) = 0. 765 END DO 766 END DO 767 768 DO i = 1, klon 769 epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. 770 epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. 771 rhow_moyen(i, 1) = rhow(i, 1) 772 END DO 773 774 DO k = 2, klev 775 DO i = 1, klon 776 epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1. 777 epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k) 778 rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* & 779 epaisseur1(i,k))/epaisseur2(i, k) 780 END DO 781 END DO 782 783 746 784 ! Choose an integration bound well above wake top 747 785 ! ----------------------------------------------------------------- 748 786 787 ! Pupper = 50000. ! melting level 788 ! Pupper = 60000. 789 ! Pupper = 80000. ! essais pour case_e 790 DO i = 1, klon 791 pupper(i) = 0.6*ph(i, 1) 792 pupper(i) = max(pupper(i), 45000.) 793 ! cc Pupper(i) = 60000. 794 END DO 795 796 749 797 ! Determine Wake top pressure (Ptop) from buoyancy integral 750 798 ! -------------------------------------------------------- 751 799 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) 800 ! -1/ Pressure of the level where dth becomes less than delta_t_min. 801 802 DO i = 1, klon 803 ptop_provis(i) = ph(i, 1) 804 END DO 805 DO k = 2, klev 806 DO i = 1, klon 807 808 ! IM v3JYG; ptop_provis(i).LT. ph(i,1) 809 810 IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. & 811 ptop_provis(i)==ph(i,1)) THEN 812 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- & 813 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 814 END IF 815 END DO 816 END DO 817 818 ! -2/ dth integral 819 820 DO i = 1, klon 821 sum_dth(i) = 0. 822 dthmin(i) = -delta_t_min 823 z(i) = 0. 824 END DO 825 826 DO k = 1, klev 827 DO i = 1, klon 828 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) 829 IF (dz(i)>0) THEN 830 z(i) = z(i) + dz(i) 831 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 832 dthmin(i) = amin1(dthmin(i), dth(i,k)) 833 END IF 834 END DO 835 END DO 836 837 ! -3/ height of triangle with area= sum_dth and base = dthmin 838 839 DO i = 1, klon 840 hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) 841 hw0(i) = amax1(hwmin, hw0(i)) 842 END DO 843 844 ! -4/ now, get Ptop 845 846 DO i = 1, klon 847 z(i) = 0. 848 ptop(i) = ph(i, 1) 849 END DO 850 851 DO k = 1, klev 852 DO i = 1, klon 853 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i)) 854 IF (dz(i)>0) THEN 855 z(i) = z(i) + dz(i) 856 ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) 857 END IF 858 END DO 859 END DO 860 861 IF (prt_level>=10) THEN 862 PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout) 863 ENDIF 864 865 866 ! -5/ Determination de ktop et kupper 867 868 DO k = klev, 1, -1 869 DO i = 1, klon 870 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 871 IF (ph(i,k+1)<pupper(i)) kupper(i) = k 872 END DO 873 END DO 874 875 ! On evite kupper = 1 et kupper = klev 876 DO i = 1, klon 877 kupper(i) = max(kupper(i), 2) 878 kupper(i) = min(kupper(i), klev-1) 879 END DO 880 881 882 ! -6/ Correct ktop and ptop 883 884 DO i = 1, klon 885 ptop_new(i) = ptop(i) 886 END DO 887 DO k = klev, 2, -1 888 DO i = 1, klon 889 IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & 890 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 891 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & 892 k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 893 END IF 894 END DO 895 END DO 896 897 DO i = 1, klon 898 ptop(i) = ptop_new(i) 899 END DO 900 901 DO k = klev, 1, -1 902 DO i = 1, klon 903 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 904 END DO 905 END DO 906 907 IF (prt_level>=10) THEN 908 PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) 763 909 ENDIF 764 910 … … 771 917 deltaqw(i, k) = 0. 772 918 d_deltatw2(i,k) = -deltatw0(i,k) 773 d_deltaqw2(i,k) = -deltaqw0(i,k)774 919 #ifdef ISO 775 920 do ixt=1,ntraciso … … 796 941 ! -------------------------------------- 797 942 798 ! Initialize sum_thv xto 1st level virt. pot. temp.943 ! Initialize sum_thvu to 1st level virt. pot. temp. 799 944 800 945 DO i = 1, klon 801 946 z(i) = 1. 802 947 dz(i) = 1. 803 sum_thv x(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)948 sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i) 804 949 sum_dth(i) = 0. 805 950 END DO … … 807 952 DO k = 1, klev 808 953 DO i = 1, klon 809 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)* RG)954 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) 810 955 IF (dz(i)>0) THEN 811 ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau812 956 z(i) = z(i) + dz(i) 813 sum_th x(i) = sum_thx(i) + thx(i, k)*dz(i)814 sum_t x(i) = sum_tx(i) + tx(i, k)*dz(i)815 sum_q x(i) = sum_qx(i) + qx(i, k)*dz(i)816 sum_thv x(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)957 sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) 958 sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) 959 sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) 960 sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i) 817 961 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 818 962 sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) 963 sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) 819 964 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) 820 965 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) … … 836 981 837 982 DO i = 1, klon 838 av_th x(i) = sum_thx(i)/hw0(i)839 av_t x(i) = sum_tx(i)/hw0(i)840 av_q x(i) = sum_qx(i)/hw0(i)841 av_thv x(i) = sum_thvx(i)/hw0(i)983 av_thu(i) = sum_thu(i)/hw0(i) 984 av_tu(i) = sum_tu(i)/hw0(i) 985 av_qu(i) = sum_qu(i)/hw0(i) 986 av_thvu(i) = sum_thvu(i)/hw0(i) 842 987 ! av_thve = sum_thve/hw0 843 988 av_dth(i) = sum_dth(i)/hw0(i) 844 989 av_dq(i) = sum_dq(i)/hw0(i) 990 av_rho(i) = sum_rho(i)/hw0(i) 845 991 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 846 992 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 847 993 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) 850 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 994 wape(i) = -rg*hw0(i)*(av_dth(i)+ & 995 epsim1*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) 996 997 END DO 855 998 856 999 ! 2.2 Prognostic variable update … … 879 1022 DO i = 1, klon 880 1023 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_targ886 ENDIF !! (wape(i)<0.)887 ENDDO888 !889 IF (iflag_wk_pop_dyn == 3) THEN890 DO i = 1, klon891 IF (wape(i)<0.) THEN892 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_targ896 ENDIF !! (wape(i)<0.)897 ENDDO898 ENDIF !! (iflag_wk_pop_dyn == 3)899 900 DO i = 1, klon901 IF (wape(i)<0.) THEN902 1024 wape(i) = 0. 903 1025 cstar(i) = 0. 904 1026 hw(i) = hwmin 1027 !jyg< 1028 !! sigmaw(i) = amax1(sigmad, sigd_con(i)) 1029 sigmaw_targ = max(sigmad, sigd_con(i)) 1030 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 1031 sigmaw(i) = sigmaw_targ 1032 !>jyg 905 1033 fip(i) = 0. 906 1034 gwake(i) = .FALSE. … … 911 1039 END IF 912 1040 END DO 913 ! 1041 914 1042 915 1043 ! Check qx and qw positivity 916 1044 ! -------------------------- 917 1045 DO i = 1, klon 918 q0_min(i) = min((q b(i,1)-sigmaw(i)*deltaqw(i,1)), &919 (q b(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))1046 q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), & 1047 (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1))) 920 1048 END DO 921 1049 DO k = 2, klev 922 1050 DO i = 1, klon 923 q1_min(i) = min((q b(i,k)-sigmaw(i)*deltaqw(i,k)), &924 (q b(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))1051 q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), & 1052 (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k))) 925 1053 IF (q1_min(i)<=q0_min(i)) THEN 926 1054 q0_min(i) = q1_min(i) … … 932 1060 ok_qx_qw(i) = q0_min(i) >= 0. 933 1061 alpha(i) = 1. 934 alpha_tot(i) = 1.935 1062 END DO 936 1063 … … 945 1072 ! ----------------- 946 1073 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 ! ------------------------------------------------------------------------ 1074 nsub = 10 1075 dtimesub = dtime/nsub 1076 1077 ! ------------------------------------------------------------ 1078 DO isubstep = 1, nsub 1079 ! ------------------------------------------------------------ 1080 959 1081 ! wk_adv is the logical flag enabling wake evolution in the time advance 960 1082 ! loop … … 965 1087 PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', & 966 1088 isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) 967 968 1089 ENDIF 969 1090 970 1091 ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement 971 ! n egatif de ktop akupper --------972 ! cc On calcule pour cela une densit ewdens0 pour laquelle on1092 ! négatif de ktop à kupper -------- 1093 ! cc On calcule pour cela une densité wdens0 pour laquelle on 973 1094 ! aurait un entrainement nul --- 974 1095 !jyg< … … 977 1098 ! des descentes unsaturees. Nous faisons alors l'hypothese que la 978 1099 ! convection profonde cree directement de nouvelles poches, sans passer 979 ! par les thermiques. La nouvelle valeur de wdens est alors impos ee.1100 ! par les thermiques. La nouvelle valeur de wdens est alors imposée. 980 1101 981 1102 DO i = 1, klon … … 983 1104 ! c $ isubstep,wk_adv(i),cstar(i),wape(i) 984 1105 IF (wk_adv(i) .AND. cstar(i)>0.01) THEN 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 omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & 1107 rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) 991 1108 wdens0 = (sigmaw(i)/(4.*3.14))* & 992 1109 ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2) … … 997 1114 IF (wdens(i)<=wdens0*1.1) THEN 998 1115 IF (iflag_wk_pop_dyn >= 1) THEN 999 d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)1000 1116 d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i) 1001 1117 ENDIF … … 1005 1121 END DO 1006 1122 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 !!-------------------------------------------------------- 1123 DO i = 1, klon 1124 IF (wk_adv(i)) THEN 1125 gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) 1126 rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i))) 1127 !jyg< 1128 !! sigmaw(i) = amin1(sigmaw(i), sigmaw_max) 1129 sigmaw_targ = min(sigmaw(i), sigmaw_max) 1130 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 1131 sigmaw(i) = sigmaw_targ 1132 !>jyg 1133 END IF 1134 END DO 1135 1136 IF (iflag_wk_pop_dyn >= 1) THEN 1137 ! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0. 1138 ! Here, it has to be set to zero. 1139 death_rate(:) = 0. 1140 1141 IF (iflag_wk_act ==2) THEN 1013 1142 DO i = 1, klon 1014 1143 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 1022 DO i = 1, klon 1023 IF (wk_adv(i)) THEN 1024 sigmaw_targ = min(sigmaw(i), sigmaw_max) 1025 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) 1026 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 1027 sigmaw(i) = sigmaw_targ 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 !!-------------------------------------------------------- 1144 wape1_act(i) = abs(cin(i)) 1145 wape2_act(i) = 2.*wape1_act(i) + 1. 1146 act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) )) 1147 ENDIF ! (wk_adv(i)) 1148 ENDDO 1149 ENDIF ! (iflag_wk_act ==2) 1150 1151 1035 1152 DO i = 1, klon 1036 1153 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 !!-------------------------------------------------------- 1043 1044 IF (iflag_wk_pop_dyn >= 1) THEN 1045 ! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0. 1046 ! Here, it has to be set to zero. 1047 death_rate(:) = 0. 1048 ENDIF 1049 1050 IF (iflag_wk_pop_dyn >= 3) THEN 1051 DO i = 1, klon 1052 IF (wk_adv(i)) THEN 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 1154 !! tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.) 1155 tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.) 1156 tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub) 1157 drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / & 1158 (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i)) 1159 !! (1 - 2*sigmaw(i)*(1.-f_shear(i))) 1160 drdt_pos=max(drdt(i),0.) 1161 1162 !! d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) & 1163 !! - wdens(i)*tau_wk_inv_min & 1164 !! - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub 1165 d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub 1166 d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min - & 1167 2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub 1168 d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i)) 1169 1170 !! d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) & 1171 !! + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) & 1172 !! - sigmaw(i)*tau_wk_inv_min )*dtimesub 1173 d_sig_gen(i) = wgen(i)*aa0 1174 d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min 1175 !! d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos 1176 d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos 1177 d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub 1178 d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) 1057 1179 ENDIF 1058 1180 ENDDO 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 1181 1182 IF (prt_level >= 10) THEN 1183 print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', & 1184 cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) 1185 print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', & 1186 wdens(1), awdens(1), act(1), d_awdens(1) 1187 print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', & 1188 wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1) 1189 print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & 1190 d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) 1191 ENDIF 1075 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 ELSE ! (iflag_wk_pop_dyn >= 1) 1194 1125 1195 ! cc nrlmd 1126 1196 1127 1197 DO i = 1, klon 1128 1198 IF (wk_adv(i)) THEN 1129 1130 ! cc nrlmd Introduction du taux de mortalite des poches et 1199 ! cc nrlmd Introduction du taux de mortalité des poches et 1131 1200 ! test sur sigmaw_max=0.4 1132 1201 ! cc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub … … 1153 1222 END DO 1154 1223 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 1224 ENDIF ! (iflag_wk_pop_dyn >= 1) 1225 1226 1167 1227 ! calcul de la difference de vitesse verticale poche - zone non perturbee 1168 1228 ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg 1169 1229 ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit 1170 ! IM 060208 au niveau k=1.. .1230 ! IM 060208 au niveau k=1..? 1171 1231 !JYG 161013 Correction : maintenant omg est dimensionne a klev. 1172 1232 DO k = 1, klev … … 1196 1256 DO i = 1, klon 1197 1257 IF (wk_adv(i) .AND. k<=ktop(i)) THEN 1198 dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)* RG)1258 dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg) 1199 1259 z(i) = z(i) + dz(i) 1200 1260 dp_deltomg(i, k) = dp_deltomg(i, 1) … … 1206 1266 DO i = 1, klon 1207 1267 IF (wk_adv(i)) THEN 1208 dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))* RG)1268 dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg) 1209 1269 ztop(i) = z(i) + dztop(i) 1210 1270 omgtop(i) = dp_deltomg(i, 1)*ztop(i) … … 1224 1284 DO i = 1, klon 1225 1285 IF (wk_adv(i)) THEN 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 omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i) 1287 dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) 1229 1288 END IF 1230 1289 END DO … … 1233 1292 DO i = 1, klon 1234 1293 IF (wk_adv(i) .AND. k<=ktop(i)) THEN 1235 omg(i, k) = -rho(i, k)* RG*omg(i, k)1294 omg(i, k) = -rho(i, k)*rg*omg(i, k) 1236 1295 dp_deltomg(i, k) = dp_deltomg(i, 1) 1237 1296 END IF … … 1243 1302 DO i = 1, klon 1244 1303 IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN 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 omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & 1305 rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) 1251 1306 dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ & 1252 1307 (ptop(i)-pupper(i)) … … 1255 1310 1256 1311 ! c DO i=1,klon 1257 ! c print*,'Pente entre 0 et kupper (r eference)'1312 ! c print*,'Pente entre 0 et kupper (référence)' 1258 1313 ! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) 1259 1314 ! c print*,'Pente entre ktop et kupper' … … 1331 1386 IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN 1332 1387 dth(i, k) = deltatw(i, k)/ppi(i, k) 1333 th1(i, k) = th b(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area1334 th2(i, k) = th b(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake1335 q1(i, k) = q b(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area1336 q2(i, k) = q b(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake1388 th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area 1389 th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake 1390 q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area 1391 q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake 1337 1392 #ifdef ISO 1338 1393 do ixt=1,ntraciso 1339 xt1(ixt,i,k) = xt b(ixt,i,k) -sigmaw(i) *deltaxtw(ixt,i,k) ! undisturbed area1340 xt2(ixt,i,k) = xt b(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake1394 xt1(ixt,i,k) = xte(ixt,i,k) -sigmaw(i) *deltaxtw(ixt,i,k) ! undisturbed area 1395 xt2(ixt,i,k) = xte(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake 1341 1396 enddo 1342 1397 #endif … … 1416 1471 END DO 1417 1472 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)) 1473 IF (prt_level>=10) THEN 1474 PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,klev) 1475 PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,klev) 1476 PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,klev) 1477 PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,klev) 1424 1478 ENDIF 1425 1479 … … 1437 1491 (1.-alpha_up(i,k))*omgbdth(i,k)- & 1438 1492 alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k) 1439 ! print*,'d_d ,k_ptop_provis(i)eltatw=', k, d_deltatw(i,k)1493 ! print*,'d_deltatw=', k, d_deltatw(i,k) 1440 1494 1441 1495 d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & … … 1471 1525 ! C 1472 1526 ! ----------------------------------------------------------------- 1473 d_t b(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &1527 d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- & 1474 1528 rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/ & 1475 1529 (ph(i,k)-ph(i,k+1)) & … … 1477 1531 (ph(i,k)-ph(i,k+1)) )*ppi(i, k) 1478 1532 1479 d_q b(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- &1533 d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- & 1480 1534 rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ & 1481 1535 (ph(i,k)-ph(i,k+1)) & … … 1484 1538 #ifdef ISO 1485 1539 do ixt=1,ntraciso 1486 d_xt b(ixt,i,k) = dtimesub*( &1540 d_xte(ixt,i,k) = dtimesub*( & 1487 1541 ( rre1(i)*omg(i,k )*sigmaw(i) *d_xt1(ixt,i,k) & 1488 1542 -rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1) ) & … … 1494 1548 #endif 1495 1549 ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN 1496 d_t b(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_q b(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))1550 d_te(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k) 1551 1552 d_qe(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1))) 1499 1553 1500 1554 #ifdef ISO 1501 1555 do ixt=1,ntraciso 1502 d_xt b(ixt,i,k) = dtimesub*( &1556 d_xte(ixt,i,k) = dtimesub*( & 1503 1557 ( rre1(i)*omg(i,k )*sigmaw(i) *d_xt1(ixt,i,k) & 1504 1558 /(Ph(i,k)-Ph(i,k+1))) & … … 1550 1604 1551 1605 1552 ! Coefficient de r epartition1606 ! Coefficient de répartition 1553 1607 1554 1608 crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ & 1555 1609 (ph(i,kupper(i))-ph(i,1)) 1556 1610 crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ & 1557 (p h(i,1)-ph(i,kupper(i)))1611 (p(i,1)-ph(i,kupper(i))) 1558 1612 1559 1613 … … 1594 1648 ! 1595 1649 1596 ! cc nrlmd Prise en compte du taux de mortalit e1597 ! cc D efinitions de entr, detr1650 ! cc nrlmd Prise en compte du taux de mortalité 1651 ! cc Définitions de entr, detr 1598 1652 !jyg< 1599 1653 !! detr(i, k) = 0. … … 1605 1659 sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) 1606 1660 !>jyg 1607 wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)1608 1609 ! cc wkspread(i,k) =1661 spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i) 1662 1663 ! cc spread(i,k) = 1610 1664 ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ 1611 1665 ! cc $ sigmaw(i) 1612 1666 1613 1667 1614 ! ajout d'un effet onde de gravit e-Tgw(k)*deltatw(k) 03/02/06 YU1668 ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU 1615 1669 ! Jingmei 1616 1670 … … 1624 1678 ! Sans GW 1625 1679 1626 ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)- wkspread(k)*deltatw(k))1680 ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) 1627 1681 1628 1682 ! GW formule 1 1629 1683 1630 1684 ! deltatw(k) = deltatw(k)+dtimesub* 1631 ! $ (ff+dtKE(k) - wkspread(k)*deltatw(k)-Tgw(k)*deltatw(k))1685 ! $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) 1632 1686 1633 1687 ! GW formule 2 … … 1689 1743 ! Scale tendencies so that water vapour remains positive in w and x. 1690 1744 1691 CALL wake_vec_modulation(klon, klev, wk_adv, epsilon _loc, qb, d_qb, deltaqw, &1745 CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, & 1692 1746 d_deltaqw, sigmaw, d_sigmaw, alpha) 1693 !1694 ! Alpha_tot = Product of all the alpha's1695 DO i = 1, klon1696 IF (wk_adv(i)) THEN1697 alpha_tot(i) = alpha_tot(i)*alpha(i)1698 END IF1699 END DO1700 1747 1701 1748 ! cc nrlmd … … 1708 1755 DO i = 1, klon 1709 1756 IF (wk_adv(i) .AND. k<=kupper(i)) THEN 1710 d_t b(i, k) = alpha(i)*d_tb(i, k)1711 d_q b(i, k) = alpha(i)*d_qb(i, k)1757 d_te(i, k) = alpha(i)*d_te(i, k) 1758 d_qe(i, k) = alpha(i)*d_qe(i, k) 1712 1759 d_deltatw(i, k) = alpha(i)*d_deltatw(i, k) 1713 1760 d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k) … … 1715 1762 #ifdef ISO 1716 1763 do ixt=1,ntraciso 1717 d_xt b(ixt,i,k)=alpha(i)*d_xtb(ixt,i,k)1764 d_xte(ixt,i,k)=alpha(i)*d_xte(ixt,i,k) 1718 1765 d_deltaxtw(ixt,i,k)=alpha(i)*d_deltaxtw(ixt,i,k) 1719 1766 enddo !do ixt=1,ntraciso … … 1734 1781 DO i = 1, klon 1735 1782 IF (wk_adv(i) .AND. k<=kupper(i)) THEN 1736 dtls(i, k) = dtls(i, k) + d_t b(i, k)1737 dqls(i, k) = dqls(i, k) + d_q b(i, k)1783 dtls(i, k) = dtls(i, k) + d_te(i, k) 1784 dqls(i, k) = dqls(i, k) + d_qe(i, k) 1738 1785 #ifdef ISO 1739 1786 do ixt=1,ntraciso 1740 dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xt b(ixt,i,k)1787 dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xte(ixt,i,k) 1741 1788 enddo !do ixt=1,ntraciso 1742 1789 #endif … … 1772 1819 DO i = 1, klon 1773 1820 IF (wk_adv(i) .AND. k<=kupper(i)) THEN 1774 t b(i, k) = tb0(i, k) + dtls(i, k)1775 q b(i, k) = qb0(i, k) + dqls(i, k)1776 th b(i, k) = tb(i, k)/ppi(i, k)1821 te(i, k) = te0(i, k) + dtls(i, k) 1822 qe(i, k) = qe0(i, k) + dqls(i, k) 1823 the(i, k) = te(i, k)/ppi(i, k) 1777 1824 deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k) 1778 1825 deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k) 1779 1826 dth(i, k) = deltatw(i, k)/ppi(i, k) 1780 ! c print*,'k,qx,qw',k,q b(i,k)-sigmaw(i)*deltaqw(i,k)1781 ! c $ ,q b(i,k)+(1-sigmaw(i))*deltaqw(i,k)1827 ! c print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k) 1828 ! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k) 1782 1829 #ifdef ISO 1783 1830 do ixt=1,ntraciso 1784 xt b(ixt,i,k) = xtb0(ixt,i,k) + dxtls(ixt,i,k)1831 xte(ixt,i,k) = xte0(ixt,i,k) + dxtls(ixt,i,k) 1785 1832 deltaxtw(ixt,i,k) = deltaxtw(ixt,i,k)+d_deltaxtw(ixt,i,k) 1786 1833 enddo !do ixt=1,ntraciso … … 1803 1850 DO k= 1,klev 1804 1851 DO i = 1,klon 1805 call iso_verif_egalite_choix(xt b(iso_eau,i,k), &1806 q b(i,k),'wake 1379',errmax,errmaxrel)1852 call iso_verif_egalite_choix(xte(iso_eau,i,k), & 1853 qe(i,k),'wake 1379',errmax,errmaxrel) 1807 1854 enddo ! DO i = 1,klon 1808 1855 enddo ! DO k= 1,klev … … 1810 1857 if (iso_hdo.gt.0) then 1811 1858 call iso_verif_aberrant_enc_vect2D( & 1812 xt b,qb, &1813 'wake 1456, xt bapres modifs',ntraciso,klon,klev)1859 xte,qe, & 1860 'wake 1456, xte apres modifs',ntraciso,klon,klev) 1814 1861 ! call iso_verif_aberrant_enc_vect2D_ns( 1815 1862 ! : deltaxtw,deltaqw, … … 1823 1870 DO i = 1, klon 1824 1871 IF (wk_adv(i)) THEN 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 1872 awdens(i) = awdens(i) + d_awdens(i) 1873 wdens(i) = wdens(i) + d_wdens(i) 1874 d_awdens2(i) = d_awdens2(i) + d_awdens(i) 1875 d_wdens2(i) = d_wdens2(i) + d_wdens(i) 1876 END IF 1877 END DO 1878 DO i = 1, klon 1879 IF (wk_adv(i)) THEN 1880 wdens_targ = max(wdens(i),wdensmin) 1881 d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i) 1882 wdens(i) = wdens_targ 1883 ! 1884 wdens_targ = min( max(awdens(i),0.), wdens(i) ) 1885 d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i) 1886 awdens(i) = wdens_targ 1887 END IF 1888 END DO 1833 1889 DO i = 1, klon 1834 1890 IF (wk_adv(i)) THEN 1835 1891 sigmaw_targ = max(sigmaw(i),sigmad) 1836 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)1837 1892 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 1838 1893 sigmaw(i) = sigmaw_targ 1839 1894 END IF 1840 1895 END DO 1841 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1842 ! Cumulatives 1896 ENDIF ! (iflag_wk_pop_dyn >= 1) 1897 !>jyg 1898 1899 1900 ! Determine Ptop from buoyancy integral 1901 ! --------------------------------------- 1902 1903 ! - 1/ Pressure of the level where dth changes sign. 1904 1905 DO i = 1, klon 1906 IF (wk_adv(i)) THEN 1907 ptop_provis(i) = ph(i, 1) 1908 END IF 1909 END DO 1910 1911 DO k = 2, klev 1912 DO i = 1, klon 1913 IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & 1914 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 1915 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 1916 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 1917 END IF 1918 END DO 1919 END DO 1920 1921 ! - 2/ dth integral 1922 1923 DO i = 1, klon 1924 IF (wk_adv(i)) THEN !!! nrlmd 1925 sum_dth(i) = 0. 1926 dthmin(i) = -delta_t_min 1927 z(i) = 0. 1928 END IF 1929 END DO 1930 1931 DO k = 1, klev 1843 1932 DO i = 1, klon 1844 1933 IF (wk_adv(i)) THEN 1845 wdens(i) = wdens(i) + d_wdens(i) 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) 1851 END IF 1852 END DO 1853 ! Bounds 1934 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) 1935 IF (dz(i)>0) THEN 1936 z(i) = z(i) + dz(i) 1937 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 1938 dthmin(i) = amin1(dthmin(i), dth(i,k)) 1939 END IF 1940 END IF 1941 END DO 1942 END DO 1943 1944 ! - 3/ height of triangle with area= sum_dth and base = dthmin 1945 1946 DO i = 1, klon 1947 IF (wk_adv(i)) THEN 1948 hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) 1949 hw(i) = amax1(hwmin, hw(i)) 1950 END IF 1951 END DO 1952 1953 ! - 4/ now, get Ptop 1954 1955 DO i = 1, klon 1956 IF (wk_adv(i)) THEN !!! nrlmd 1957 ktop(i) = 0 1958 z(i) = 0. 1959 END IF 1960 END DO 1961 1962 DO k = 1, klev 1854 1963 DO i = 1, klon 1855 1964 IF (wk_adv(i)) THEN 1856 wdens_targ = max(wdens(i),wdensmin) 1857 d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i) 1858 d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i) 1859 wdens(i) = wdens_targ 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 1873 wdens_targ = min( max(awdens(i),0.), wdens(i) ) 1874 d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) 1875 d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i) 1876 awdens(i) = wdens_targ 1877 END IF 1878 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) 1962 ENDIF ! (iflag_wk_pop_dyn >= 1) 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 1965 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i)) 1966 IF (dz(i)>0) THEN 1967 z(i) = z(i) + dz(i) 1968 ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) 1969 ktop(i) = k 1970 END IF 1971 END IF 1972 END DO 1973 END DO 1974 1975 ! 4.5/Correct ktop and ptop 1976 1977 DO i = 1, klon 1978 IF (wk_adv(i)) THEN 1979 ptop_new(i) = ptop(i) 1980 END IF 1981 END DO 1982 1983 DO k = klev, 2, -1 1984 DO i = 1, klon 1985 ! IM v3JYG; IF (k .GE. ktop(i) 1986 IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & 1987 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 1988 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 1989 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 1990 END IF 1991 END DO 1992 END DO 1993 1994 1995 DO i = 1, klon 1996 IF (wk_adv(i)) THEN 1997 ptop(i) = ptop_new(i) 1998 END IF 1999 END DO 2000 2001 DO k = klev, 1, -1 2002 DO i = 1, klon 2003 IF (wk_adv(i)) THEN !!! nrlmd 2004 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 2005 END IF 2006 END DO 2007 END DO 1970 2008 1971 2009 ! 5/ Set deltatw & deltaqw to 0 above kupper … … 1992 2030 DO i = 1, klon 1993 2031 IF (wk_adv(i)) THEN !!! nrlmd 1994 sum_th x(i) = 0.1995 sum_t x(i) = 0.1996 sum_q x(i) = 0.1997 sum_thv x(i) = 0.2032 sum_thu(i) = 0. 2033 sum_tu(i) = 0. 2034 sum_qu(i) = 0. 2035 sum_thvu(i) = 0. 1998 2036 sum_dth(i) = 0. 1999 2037 sum_dq(i) = 0. 2038 sum_rho(i) = 0. 2000 2039 sum_dtdwn(i) = 0. 2001 2040 sum_dqdwn(i) = 0. 2002 2041 2003 av_th x(i) = 0.2004 av_t x(i) = 0.2005 av_q x(i) = 0.2006 av_thv x(i) = 0.2042 av_thu(i) = 0. 2043 av_tu(i) = 0. 2044 av_qu(i) = 0. 2045 av_thvu(i) = 0. 2007 2046 av_dth(i) = 0. 2008 2047 av_dq(i) = 0. 2048 av_rho(i) = 0. 2009 2049 av_dtdwn(i) = 0. 2010 2050 av_dqdwn(i) = 0. … … 2015 2055 ! -------------------------------------- 2016 2056 2017 ! Initialize sum_thv xto 1st level virt. pot. temp.2057 ! Initialize sum_thvu to 1st level virt. pot. temp. 2018 2058 2019 2059 DO i = 1, klon … … 2021 2061 z(i) = 1. 2022 2062 dz(i) = 1. 2023 sum_thv x(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)2063 sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i) 2024 2064 sum_dth(i) = 0. 2025 2065 END IF … … 2029 2069 DO i = 1, klon 2030 2070 IF (wk_adv(i)) THEN !!! nrlmd 2031 dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)* RG)2071 dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) 2032 2072 IF (dz(i)>0) THEN 2033 2073 z(i) = z(i) + dz(i) 2034 sum_th x(i) = sum_thx(i) + thx(i, k)*dz(i)2035 sum_t x(i) = sum_tx(i) + tx(i, k)*dz(i)2036 sum_q x(i) = sum_qx(i) + qx(i, k)*dz(i)2037 sum_thv x(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)2074 sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) 2075 sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) 2076 sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) 2077 sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i) 2038 2078 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 2039 2079 sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) 2080 sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) 2040 2081 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) 2041 2082 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) … … 2061 2102 DO i = 1, klon 2062 2103 IF (wk_adv(i)) THEN !!! nrlmd 2063 av_th x(i) = sum_thx(i)/hw0(i)2064 av_t x(i) = sum_tx(i)/hw0(i)2065 av_q x(i) = sum_qx(i)/hw0(i)2066 av_thv x(i) = sum_thvx(i)/hw0(i)2104 av_thu(i) = sum_thu(i)/hw0(i) 2105 av_tu(i) = sum_tu(i)/hw0(i) 2106 av_qu(i) = sum_qu(i)/hw0(i) 2107 av_thvu(i) = sum_thvu(i)/hw0(i) 2067 2108 av_dth(i) = sum_dth(i)/hw0(i) 2068 2109 av_dq(i) = sum_dq(i)/hw0(i) 2110 av_rho(i) = sum_rho(i)/hw0(i) 2069 2111 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 2070 2112 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 2071 2113 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 2114 wape(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + & 2115 av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) 2116 END IF 2117 END DO 2077 2118 2078 2119 ! Filter out bad wakes … … 2107 2148 !! sigmaw(i) = max(sigmad, sigd_con(i)) 2108 2149 sigmaw_targ = max(sigmad, sigd_con(i)) 2109 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)2110 2150 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 2111 2151 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_targ2116 2152 !>jyg 2117 2153 fip(i) = 0. … … 2123 2159 END IF 2124 2160 END DO 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 2161 2162 END DO ! end sub-timestep loop 2163 2138 2164 IF (prt_level>=10) THEN 2139 2165 PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', & … … 2155 2181 ! cc 2156 2182 z(i) = 0. 2157 sum_th x(i) = 0.2158 sum_t x(i) = 0.2159 sum_q x(i) = 0.2160 sum_thv x(i) = 0.2183 sum_thu(i) = 0. 2184 sum_tu(i) = 0. 2185 sum_qu(i) = 0. 2186 sum_thvu(i) = 0. 2161 2187 sum_dth(i) = 0. 2162 2188 sum_half_dth(i) = 0. 2163 2189 sum_dq(i) = 0. 2190 sum_rho(i) = 0. 2164 2191 sum_dtdwn(i) = 0. 2165 2192 sum_dqdwn(i) = 0. 2166 2193 2167 av_th x(i) = 0.2168 av_t x(i) = 0.2169 av_q x(i) = 0.2170 av_thv x(i) = 0.2194 av_thu(i) = 0. 2195 av_tu(i) = 0. 2196 av_qu(i) = 0. 2197 av_thvu(i) = 0. 2171 2198 av_dth(i) = 0. 2172 2199 av_dq(i) = 0. 2200 av_rho(i) = 0. 2173 2201 av_dtdwn(i) = 0. 2174 2202 av_dqdwn(i) = 0. … … 2185 2213 IF (ok_qx_qw(i)) THEN 2186 2214 ! cc 2187 rho(i, k) = p(i, k)/( RD*tb(i,k))2215 rho(i, k) = p(i, k)/(rd*te(i,k)) 2188 2216 IF (k==1) THEN 2189 rhoh(i, k) = ph(i, k)/( RD*tb(i,k))2217 rhoh(i, k) = ph(i, k)/(rd*te(i,k)) 2190 2218 zhh(i, k) = 0 2191 2219 ELSE 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) 2220 rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) 2221 zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) 2222 END IF 2223 the(i, k) = te(i, k)/ppi(i, k) 2224 thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 2225 tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) 2226 qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) 2227 rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) 2199 2228 dth(i, k) = deltatw(i, k)/ppi(i, k) 2200 2229 #ifdef ISO 2201 2230 do ixt=1,ntraciso 2202 xt x(ixt,i,k) = xtb(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)2231 xtu(ixt,i,k) = xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i) 2203 2232 enddo !do ixt=1,ntraciso 2204 2233 #endif … … 2211 2240 if (iso_hdo.gt.0) then 2212 2241 call iso_verif_aberrant_enc_vect2D( & 2213 xt x,qx, &2242 xtu,qu, & 2214 2243 'wake 1834, apres modifs',ntraciso,klon,klev) 2215 2244 endif … … 2220 2249 ! ----------------------------------------------------------- 2221 2250 2222 ! Initialize sum_thv xto 1st level virt. pot. temp.2251 ! Initialize sum_thvu to 1st level virt. pot. temp. 2223 2252 2224 2253 DO i = 1, klon … … 2229 2258 dz(i) = 1. 2230 2259 dz_half(i) = 1. 2231 sum_thv x(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)2260 sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i) 2232 2261 sum_dth(i) = 0. 2233 2262 END IF … … 2239 2268 IF (ok_qx_qw(i)) THEN 2240 2269 ! cc 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 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) 2271 dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*rg) 2243 2272 IF (dz(i)>0) THEN 2244 2273 z(i) = z(i) + dz(i) 2245 sum_th x(i) = sum_thx(i) + thx(i, k)*dz(i)2246 sum_t x(i) = sum_tx(i) + tx(i, k)*dz(i)2247 sum_q x(i) = sum_qx(i) + qx(i, k)*dz(i)2248 sum_thv x(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)2274 sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) 2275 sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) 2276 sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) 2277 sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i) 2249 2278 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 2250 2279 sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) 2280 sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) 2251 2281 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) 2252 2282 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) … … 2278 2308 IF (ok_qx_qw(i)) THEN 2279 2309 ! cc 2280 av_th x(i) = sum_thx(i)/hw0(i)2281 av_t x(i) = sum_tx(i)/hw0(i)2282 av_q x(i) = sum_qx(i)/hw0(i)2283 av_thv x(i) = sum_thvx(i)/hw0(i)2310 av_thu(i) = sum_thu(i)/hw0(i) 2311 av_tu(i) = sum_tu(i)/hw0(i) 2312 av_qu(i) = sum_qu(i)/hw0(i) 2313 av_thvu(i) = sum_thvu(i)/hw0(i) 2284 2314 av_dth(i) = sum_dth(i)/hw0(i) 2285 2315 av_dq(i) = sum_dq(i)/hw0(i) 2316 av_rho(i) = sum_rho(i)/hw0(i) 2286 2317 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 2287 2318 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 2288 2319 2289 wape2(i) = - RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + &2290 av_dth(i)*av_q x(i)+av_dth(i)*av_dq(i)))/av_thvx(i)2320 wape2(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + & 2321 av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) 2291 2322 END IF 2292 2323 END DO 2293 #ifdef IOPHYS_WK 2294 IF (.not.phys_sub) CALL iophys_ecrit('wape2_a',1,'wape2_a','J/kg',wape2) 2295 #endif 2324 2296 2325 2297 2326 … … 2323 2352 END DO 2324 2353 END IF 2325 #ifdef IOPHYS_WK2326 IF (.not.phys_sub) CALL iophys_ecrit('wape2_b',1,'wape2_b','J/kg',wape2)2327 #endif2328 2354 2329 2355 … … 2360 2386 !! sigmaw(i) = amax1(sigmad, sigd_con(i)) 2361 2387 sigmaw_targ = max(sigmad, sigd_con(i)) 2362 d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)2363 2388 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 2364 2389 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_targ2369 2390 !>jyg 2370 2391 fip(i) = 0. … … 2375 2396 gwake(i) = .TRUE. 2376 2397 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)) 2398 END IF 2381 2399 END DO 2382 2400 … … 2410 2428 END IF 2411 2429 END DO 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) 2430 2455 2431 ! Limitation de sigmaw 2456 2432 … … 2467 2443 DO i = 1, klon 2468 2444 kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & 2469 .NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold) 2470 !! .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin) 2445 .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin) 2471 2446 ENDDO 2472 2447 ELSE ! (iflag_wk_pop_dyn >= 1) … … 2508 2483 wape(i) = 0. 2509 2484 cstar(i) = 0. 2510 !!jyg Outside subroutine "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes2485 !!jyg Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes 2511 2486 !! hw(i) = hwmin !jyg 2512 2487 !! sigmaw(i) = sigmad !jyg 2513 2488 hw(i) = 0. !jyg 2514 2489 fip(i) = 0. 2515 !2516 2490 !! sigmaw(i) = 0. !jyg 2517 2491 sigmaw_targ = 0. 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 2492 d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) 2521 2493 sigmaw(i) = sigmaw_targ 2522 !2523 IF (iflag_wk_pop_dyn >= 3) THEN2524 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 202201242528 asigmaw(i) = sigmaw_targ2529 ELSE2530 asigmaw(i) = 0.2531 ENDIF ! (iflag_wk_pop_dyn >= 3)2532 !2533 2494 IF (iflag_wk_pop_dyn >= 1) THEN 2534 2495 !! awdens(i) = 0. 2535 2496 !! wdens(i) = 0. 2536 2497 wdens_targ = 0. 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 2498 d_wdens2(i) = wdens_targ - wdens(i) 2540 2499 wdens(i) = wdens_targ 2541 2500 wdens_targ = 0. 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 2501 d_awdens2(i) = wdens_targ - awdens(i) 2548 2502 awdens(i) = wdens_targ 2549 !! IF (iflag_wk_pop_dyn == 2) THEN2550 !! d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)2551 !! ENDIF ! (iflag_wk_pop_dyn == 2)2552 2503 ENDIF ! (iflag_wk_pop_dyn >= 1) 2553 2504 ELSE ! (kill_wake(i)) … … 2563 2514 wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout) 2564 2515 ENDIF 2565 #ifdef IOPHYS_WK2566 IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape)2567 #endif2568 2516 2569 2517 … … 2597 2545 END DO 2598 2546 !jyg< 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 2547 DO i = 1, klon 2548 d_sigmaw2(i) = d_sigmaw2(i)/dtime 2549 d_awdens2(i) = d_awdens2(i)/dtime 2550 d_wdens2(i) = d_wdens2(i)/dtime 2551 ENDDO 2640 2552 !>jyg 2641 2553 2642 RETURN 2554 2555 2556 RETURN 2643 2557 END SUBROUTINE wake 2644 2558 2645 SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon _loc, qb, d_qb, deltaqw, &2559 SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, & 2646 2560 d_deltaqw, sigmaw, d_sigmaw, alpha) 2647 2561 ! ------------------------------------------------------ 2648 ! D termination du coefficient alpha tel que les tendances2562 ! D\'etermination du coefficient alpha tel que les tendances 2649 2563 ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent 2650 2564 ! a une humidite positive dans la zone (x) et dans la zone (w). … … 2653 2567 2654 2568 ! Input 2655 REAL q b(nlon, nl), d_qb(nlon, nl)2569 REAL qe(nlon, nl), d_qe(nlon, nl) 2656 2570 REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl) 2657 2571 REAL sigmaw(nlon), d_sigmaw(nlon) … … 2664 2578 REAL alpha1(nlon) 2665 2579 REAL x, a, b, c, discrim 2666 REAL epsilon_loc 2580 REAL epsilon 2581 ! DATA epsilon/1.e-15/ 2667 2582 INTEGER i,k 2668 2583 … … 2679 2594 DO i = 1, nlon 2680 2595 IF (wk_adv(i)) THEN 2681 x = q b(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qb(i, k) + &2596 x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + & 2682 2597 (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * & 2683 2598 (deltaqw(i,k)+d_deltaqw(i,k)) 2684 2599 a = -d_sigmaw(i)*d_deltaqw(i, k) 2685 b = d_q b(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &2600 b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - & 2686 2601 deltaqw(i, k)*d_sigmaw(i) 2687 c = q b(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc2602 c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon 2688 2603 discrim = b*b - 4.*a*c 2689 2604 ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim 2690 IF (a+b>=0.) THEN !! Condition suffisante pour la positivit ede ovap2605 IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap 2691 2606 alpha1(i) = 1. 2692 2607 ELSE … … 2714 2629 END SUBROUTINE wake_vec_modulation 2715 2630 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_pupper2723 USE lmdz_wake_ini , ONLY : RG2724 USE lmdz_wake_ini , ONLY : hwmin2725 USE lmdz_wake_ini , ONLY : iflag_wk_new_ptop, wk_delta_t_min, wk_frac_int_delta_t2726 USE lmdz_wake_ini , ONLY : wk_int_delta_t_min2727 2728 IMPLICIT NONE2729 2730 INTEGER, INTENT(IN) :: klon,klev2731 REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: ph, p2732 REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: rho2733 LOGICAL, DIMENSION (klon) , INTENT(IN) :: wk_adv2734 REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: dth2735 REAL, INTENT(IN) :: delta_t_min_in2736 2737 2738 REAL, DIMENSION (klon) , INTENT(OUT) :: hw_2739 REAL, DIMENSION (klon) , INTENT(OUT) :: ptop2740 INTEGER, DIMENSION (klon) , INTENT(OUT) :: Ktop2741 REAL, DIMENSION (klon) , INTENT(OUT) :: pupper2742 INTEGER, DIMENSION (klon) , INTENT(OUT) :: kupper2743 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,k2748 2749 LOGICAL, DIMENSION (klon) :: wk_active2750 REAL :: delta_t_min2751 REAL, DIMENSION (klon) :: dthmin2752 REAL, DIMENSION (klon) :: ptop_provis,ptop_new2753 REAL, DIMENSION (klon) :: z, dz2754 REAL, DIMENSION (klon) :: sum_dth2755 2756 INTEGER, DIMENSION (klon) :: k_ptop_provis2757 REAL, DIMENSION (klon) :: zk_ptop_provis2758 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=02767 2768 2769 2770 !INTEGER, SAVE :: compte=02771 2772 ! LJYF : a priori z, dz sum_dth sont aussi des variables internes2773 ! Les eliminer apres verification convergence numerique2774 2775 !compte=compte+12776 !print*,'compte=',compte2777 2778 ! Determine Ptop from buoyancy integral2779 ! ---------------------------------------2780 2781 ! - 1/ Pressure of the level where dth changes sign.2782 !print*,'WAKE LJYF'2783 2784 2785 if (iflag_wk_new_ptop==0) then2786 delta_t_min=delta_t_min_in2787 else2788 delta_t_min=wk_delta_t_min2789 endif2790 2791 DO i = 1, klon2792 ptop_provis(i) = ph(i, 1)2793 k_ptop_provis(i) = 12794 END DO2795 2796 DO k = 2, klev2797 DO i = 1, klon2798 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) THEN2800 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN2801 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) = k2804 END IF2805 END DO2806 END DO2807 2808 2809 2810 ! - 2/ dth integral2811 2812 DO i = 1, klon2813 IF (wk_adv(i)) THEN !!! nrlmd2814 sum_dth(i) = 0.2815 dthmin(i) = -delta_t_min2816 z(i) = 0.2817 END IF2818 END DO2819 2820 DO k = 1, klev2821 DO i = 1, klon2822 IF (wk_adv(i)) THEN2823 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG)2824 IF (dz(i)>0) THEN2825 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 IF2829 END IF2830 END DO2831 END DO2832 2833 ! - 3/ height of triangle with area= sum_dth and base = dthmin2834 2835 DO i = 1, klon2836 IF (wk_adv(i)) THEN2837 hw_(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)2838 hw_(i) = amax1(hwmin, hw_(i))2839 END IF2840 END DO2841 2842 ! - 4/ now, get Ptop2843 2844 DO i = 1, klon2845 IF (wk_adv(i)) THEN !!! nrlmd2846 ktop(i) = 02847 z(i) = 0.2848 END IF2849 END DO2850 2851 DO k = 1, klev2852 DO i = 1, klon2853 IF (wk_adv(i)) THEN2854 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw_(i)-z(i))2855 IF (dz(i)>0) THEN2856 z(i) = z(i) + dz(i)2857 ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i)2858 ktop(i) = k2859 END IF2860 END IF2861 END DO2862 END DO2863 2864 ! 4.5/Correct ktop and ptop2865 2866 DO i = 1, klon2867 ptop_new(i) = ptop(i)2868 END DO2869 2870 DO k = klev, 2, -12871 DO i = 1, klon2872 ! 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) THEN2875 dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN2876 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 IF2879 END DO2880 END DO2881 2882 2883 DO i = 1, klon2884 ptop(i) = ptop_new(i)2885 END DO2886 2887 DO k = klev, 1, -12888 DO i = 1, klon2889 IF (wk_adv(i)) THEN !!! nrlmd2890 IF (ph(i,k+1)<ptop(i)) ktop(i) = k2891 END IF2892 END DO2893 END DO2894 2895 ! IF (prt_level>=10) THEN2896 ! PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)2897 ! ENDIF2898 2899 ! -----------------------------------------------------------------------2900 ! nouveau calcul de hw et ptop2901 ! -----------------------------------------------------------------------2902 !if (iflag_wk_new_ptop>0) then2903 do i=1,klon2904 ptop1(i)=ph(i,1)2905 ktop1(i)=12906 h_zzz(i)=0.2907 enddo2908 2909 IF (iflag_wk_new_ptop/=0) THEN2910 2911 int_dth(1:klon,1:klev+1)=0.2912 DO i = 1, klon2913 IF (wk_adv(i)) THEN2914 int_dth(i,1) = 0.2915 END IF2916 END DO2917 2918 if (abs(iflag_wk_new_ptop) == 1 ) then2919 DO k = 2, klev+12920 Do i = 1, klon2921 IF (wk_adv(i)) THEN2922 if (k<=k_ptop_provis(i)) then2923 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 else2926 ddd=0.2927 endif2928 int_dth(i,k) = int_dth(i,k-1) + ddd2929 !ELSE2930 ! int_dth(i,k) = 0.2931 END IF2932 END DO2933 END DO2934 else2935 k_ptop_provis(:)=klev+12936 dthmin(:)=dth(:,1)2937 ! calcul de l'int??grale de dT * dP jusqu'au dernier2938 ! niveau avec dT<0. (en s'assurant qu'on a bien un2939 ! dT negatif plus bas)2940 DO k = 1, klev2941 DO i = 1, klon2942 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.) then2945 if (k>=k_ptop_provis(i)) then2946 ddd=0.2947 else if (dth(i,k)>=0.) then2948 ddd=0.2949 k_ptop_provis(i)=k+12950 endif2951 endif2952 int_dth(i,k+1) = int_dth(i,k)+ ddd2953 ENDDO2954 ENDDO2955 2956 DO i = 1, klon2957 if ( k_ptop_provis(i)==klev+1 .or. .not. wk_adv(i)) then2958 k_ptop_provis(i)=12959 endif2960 ENDDO2961 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_min2968 do i=1,klon2969 if (int_dth(i,k_ptop_provis(i)) > -wk_int_delta_t_min .or. k_ptop_provis(i)==1) then2970 !if (1==0) then2971 wk_active(i)=.false.2972 ptop(i)=ph(i,1)2973 ktop(i)=12974 hw_(i)=0.2975 else2976 wk_active(i)=wk_adv(i)2977 endif2978 enddo2979 2980 DO i=1,klon2981 IF (wk_active(i)) THEN2982 frac_int_dth(i)=wk_frac_int_delta_t*int_dth(i,k_ptop_provis(i))2983 ENDIF2984 ENDDO2985 DO k = 1,klev2986 DO i =1, klon2987 ! print*,ipas,'yyy ',k,int_dth(i,k),frac_int_dth(i)2988 IF (wk_active(i)) THEN2989 IF (int_dth(i,k)>=frac_int_dth(i)) THEN2990 ktop1(i) = min(k, k_ptop_provis(i))2991 !ktop1(i) = k2992 !print*,ipas,'yyy ktop1= ',ktop12993 ENDIF2994 ENDIF2995 END DO2996 END DO2997 !print*, 'LAMINE'2998 2999 DO i = 1, klon3000 IF (wk_active(i)) THEN3001 !print*, ipas,'xxx1, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ',ktop13002 ddd=int_dth(i,ktop1(i)+1)-int_dth(i,ktop1(i))3003 if (ddd==0.) then3004 omega(i)=0.3005 else3006 omega(i) = (frac_int_dth(i) - int_dth(i,ktop1(i)))/ddd3007 endif3008 !! print*,'OMEGA ',omega(i)3009 END IF3010 END DO3011 3012 !! print*, 'xxx'3013 DO i = 1, klon3014 IF (wk_active(i)) THEN3015 ! 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 IF3021 END DO3022 3023 DO i=1, klon3024 IF (wk_active(i)) THEN3025 zzz(i, 1) = 03026 END IF3027 END DO3028 DO k = 1, klev3029 DO i = 1, klon3030 IF (wk_active(i)) THEN3031 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 IF3034 END DO3035 END DO3036 3037 DO i =1, klon3038 IF (wk_active(i)) THEN3039 h_zzz(i) = max((1- omega(i))*zzz(i,ktop1(i)) + omega(i)*zzz(i,ktop1(i)+1), hwmin)3040 END IF3041 END DO3042 3043 3044 ENDIF ! (iflag_wk_new_ptop/=0)3045 3046 !if (iflag_wk_new_ptop==2) then3047 IF (iflag_wk_new_ptop>0) THEN3048 do i=1,klon3049 ptop(i)=ptop1(i)3050 ktop(i)=ktop1(i)3051 hw_(i)=h_zzz(i)3052 enddo3053 3054 !endif3055 ENDIF3056 3057 kupper = 03058 3059 IF (wk_pupper<1.) THEN3060 ! Choose an integration bound well above wake top3061 ! -----------------------------------------------------------------3062 3063 ! Pupper = 50000. ! melting level3064 ! Pupper = 60000.3065 ! Pupper = 80000. ! essais pour case_e3066 DO i = 1, klon3067 ! 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 DO3072 3073 ELSE3074 DO i=1, klon3075 ! 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 DO3079 END IF3080 3081 ! -5/ Determination de kupper3082 3083 DO k = klev, 1, -13084 DO i = 1, klon3085 IF (ph(i,k+1)<pupper(i)) kupper(i) = k3086 END DO3087 END DO3088 3089 ! On evite kupper = 1 et kupper = klev3090 DO i = 1, klon3091 kupper(i) = max(kupper(i), 2)3092 kupper(i) = min(kupper(i), klev-1)3093 END DO3094 !---------- FIN nouveau calcul hw et ptop -------------------------------------3095 3096 IF (iflag_wk_new_ptop==999) THEN3097 DO i = 1, klon3098 hw_(i)=0.3099 ptop(i)=ph(i,1)3100 Ktop(i)=13101 pupper(i)=ph(i,2)3102 kupper(i)=23103 h_zzz(i)=0.3104 Ptop1(i)=ph(i,1)3105 ENDDO3106 ENDIF3107 3108 zk_ptop_provis=k_ptop_provis3109 3110 RETURN3111 END SUBROUTINE pkupper3112 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_ini3126 USE lmdz_wake_ini , ONLY : prt_level,RG3127 USE lmdz_wake_ini , ONLY : stark, wdens_ref3128 USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa03129 !! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin3130 USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn3131 USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max3132 3133 IMPLICIT NONE3134 3135 INTEGER, INTENT(IN) :: klon,klev3136 LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv3137 REAL, INTENT(IN) :: dtime3138 REAL, INTENT(IN) :: dtimesub3139 REAL, INTENT(IN) :: wdensmin3140 REAL, DIMENSION (klon), INTENT(IN) :: wgen3141 REAL, DIMENSION (klon), INTENT(IN) :: wdens3142 REAL, DIMENSION (klon), INTENT(IN) :: awdens3143 REAL, DIMENSION (klon), INTENT(IN) :: sigmaw3144 REAL, DIMENSION (klon), INTENT(IN) :: cstar3145 REAL, DIMENSION (klon), INTENT(IN) :: cin, wape3146 REAL, DIMENSION (klon), INTENT(IN) :: f_shear3147 INTEGER, INTENT(IN) :: iflag_wk_act3148 3149 3150 !3151 3152 ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields3153 ! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)3154 REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk3155 REAL, DIMENSION (klon), INTENT(OUT) :: gfl3156 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_awdens, d_wdens3157 REAL, DIMENSION (klon), INTENT(OUT) :: drdt3158 ! Some components of the tendencies of state variables3159 REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd3160 REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_spread3161 REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd3162 REAL, INTENT(OUT) :: d_wdens_targ, d_sigmaw_targ3163 3164 3165 REAL :: delta_t_min3166 INTEGER :: i, k3167 REAL :: wdens03168 ! IM 0802083169 LOGICAL, DIMENSION (klon) :: gwake3170 3171 ! Variables liees a la dynamique de population3172 REAL, DIMENSION(klon) :: act3173 REAL, DIMENSION(klon) :: tau_wk_inv3174 REAL, DIMENSION(klon) :: wape1_act, wape2_act3175 LOGICAL, DIMENSION (klon) :: kill_wake3176 REAL :: drdt_pos3177 REAL :: tau_wk_inv_min3178 3179 3180 3181 IF (iflag_wk_act == 0) THEN3182 act(:) = 0.3183 ELSEIF (iflag_wk_act == 1) THEN3184 act(:) = 1.3185 ELSEIF (iflag_wk_act ==2) THEN3186 DO i = 1, klon3187 IF (wk_adv(i)) THEN3188 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 ENDDO3193 ENDIF ! (iflag_wk_act ==2)3194 3195 DO i = 1, klon3196 IF (wk_adv(i)) THEN3197 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 IF3200 END DO3201 3202 DO i = 1, klon3203 IF (wk_adv(i)) THEN3204 !! 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) )*dtimesub3215 !jyg+mlt<3216 d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub3217 d_dens_gen(i) = wgen(i)3218 d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min3219 d_dens_col(i) = -2.*wdens(i)*gfl(i)*drdt_pos3220 d_dens_gen(i) = d_dens_gen(i)*dtimesub3221 d_dens_death(i) = d_dens_death(i)*dtimesub3222 d_dens_col(i) = d_dens_col(i)*dtimesub3223 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 )*dtimesub3227 !>jyg+mlt3228 !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_targ3234 !! d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))3235 !>jyg3236 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 )*dtimesub3241 d_sig_gen(i) = wgen(i)*aa03242 d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min3243 !!3244 3245 d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos3246 d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos3247 d_sig_spread(i) = gfl(i)*cstar(i)3248 d_sig_gen(i) = d_sig_gen(i)*dtimesub3249 d_sig_death(i) = d_sig_death(i)*dtimesub3250 d_sig_col(i) = d_sig_col(i)*dtimesub3251 d_sig_spread(i) = d_sig_spread(i)*dtimesub3252 d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)3253 !>jyg+mlt3254 !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_targ3261 !! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))3262 !>jyg3263 ENDIF3264 ENDDO3265 3266 IF (prt_level >= 10) THEN3267 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 ENDIF3276 3277 RETURN3278 END SUBROUTINE wake_popdyn_13279 3280 SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &3281 wdensmin, &3282 sigmaw, wdens, awdens, & !! states variables3283 gfl, cstar, cin, wape, rad_wk, &3284 d_sigmaw, d_wdens, d_awdens, & !! tendences3285 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_ini3293 USE lmdz_wake_ini , ONLY : prt_level,RG3294 USE lmdz_wake_ini , ONLY : stark, wdens_ref3295 USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa03296 !! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin3297 USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn3298 USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max3299 3300 IMPLICIT NONE3301 3302 INTEGER, INTENT(IN) :: klon,klev3303 LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv3304 REAL, INTENT(IN) :: dtimesub3305 REAL, INTENT(IN) :: wdensmin3306 REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes3307 REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes3308 REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area3309 REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area3310 REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes3311 REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre3312 3313 !3314 REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = wake radius3315 REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front lenght per unit area3316 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_wdens, d_awdens3317 REAL, DIMENSION (klon), INTENT(OUT) :: cont_fact !! RM facteur de contact = 2 pi * rad * C*3318 ! Some components of the tendencies of state variables3319 REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd3320 REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd3321 REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd3322 3323 3324 !! internal variables3325 3326 INTEGER :: i, k3327 REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes3328 REAL :: tau_wk_inv_min3329 REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes3330 REAL :: d_wdens_targ, d_sigmaw_targ3331 3332 3333 !! Equations3334 !! dD/dt = B - (D-A)/tau - f D^23335 !! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^23336 !! 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, klon3342 IF (wk_adv(i)) THEN3343 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 IF3346 END DO3347 3348 3349 DO i = 1, klon3350 IF (wk_adv(i)) THEN3351 !! 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_cv3355 !! 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) ! bug3358 !! 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)*aa03362 d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min3363 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)*dtimesub3367 d_sig_death(i) = d_sig_death(i)*dtimesub3368 d_sig_col(i) = d_sig_col(i)*dtimesub3369 d_sig_spread(i) = d_sig_spread(i)*dtimesub3370 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_targ3378 !! 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_min3383 d_dens_col(i) = - cont_fact(i)*wdens(i)**23384 !3385 d_dens_gen(i) = d_dens_gen(i)*dtimesub3386 d_dens_death(i) = d_dens_death(i)*dtimesub3387 d_dens_col(i) = d_dens_col(i)*dtimesub3388 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))**23393 d_adens_acol(i) = - cont_fact(i)*awdens(i)**23394 !3395 d_adens_death(i) = d_adens_death(i)*dtimesub3396 d_adens_icol(i) = d_adens_icol(i)*dtimesub3397 d_adens_acol(i) = d_adens_acol(i)*dtimesub3398 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_targ3405 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_targ3410 3411 3412 3413 ENDIF3414 ENDDO3415 3416 IF (prt_level >= 10) THEN3417 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 ENDIF3424 sigmaw=sigmaw+d_sigmaw3425 wdens=wdens+d_wdens3426 awdens=awdens+d_awdens3427 3428 RETURN3429 END SUBROUTINE wake_popdyn_23430 3431 SUBROUTINE wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &3432 wdensmin, &3433 sigmaw, asigmaw, wdens, awdens, & !! state variables3434 gfl, agfl, cstar, cin, wape, &3435 rad_wk, arad_wk, irad_wk, &3436 d_sigmaw, d_asigmaw, d_wdens, d_awdens, & !! tendencies3437 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_ini3445 USE lmdz_wake_ini , ONLY : prt_level,RG3446 USE lmdz_wake_ini , ONLY : stark, wdens_ref3447 USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa03448 !! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin3449 USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn3450 USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max3451 USE lmdz_wake_ini , ONLY : smallestreal3452 3453 IMPLICIT NONE3454 3455 INTEGER, INTENT(IN) :: klon,klev3456 LOGICAL, INTENT(IN) :: phys_sub3457 LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv3458 REAL, INTENT(IN) :: dtimesub3459 REAL, INTENT(IN) :: wdensmin3460 REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes3461 REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes3462 REAL, DIMENSION (klon), INTENT(INOUT) :: asigmaw !! sigma = fractional area of active wakes3463 REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area3464 REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area3465 REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes3466 REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre3467 3468 !3469 REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = mean wake radius3470 REAL, DIMENSION (klon), INTENT(OUT) :: arad_wk !! r_A = wake radius of active wakes3471 REAL, DIMENSION (klon), INTENT(OUT) :: irad_wk !! r_I = wake radius of inactive wakes3472 REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front length per unit area3473 REAL, DIMENSION (klon), INTENT(OUT) :: agfl !! LgA = gust front length of active wakes3474 !! per unit area3475 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_asigmaw, d_wdens, d_awdens3476 ! Some components of the tendencies of state variables3477 REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd3478 REAL, DIMENSION (klon), INTENT(OUT) :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd3479 REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd3480 REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_acol, d_adens_icol, d_adens_bnd3481 3482 3483 !! internal variables3484 3485 INTEGER :: i, k3486 REAL, DIMENSION (klon) :: iwdens, isigmaw !! inactive wake density and fractional area3487 !! REAL, DIMENSION (klon) :: d_arad, d_irad3488 REAL, DIMENSION (klon) :: igfl !! LgI = gust front length of inactive wakes3489 !! per unit area3490 REAL, DIMENSION (klon) :: s_wk !! mean area of individual wakes3491 REAL, DIMENSION (klon) :: as_wk !! mean area of individual active wakes3492 REAL, DIMENSION (klon) :: is_wk !! mean area of individual inactive wakes3493 REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes3494 REAL :: tau_wk_inv_min3495 REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes3496 REAL :: d_wdens_targ, d_sigmaw_targ3497 3498 3499 !! Equations3500 !! ---------3501 !! Gust fronts:3502 !! Lg_A = 2 pi r_A A3503 !! Lg_I = 2 pi r_I I3504 !! Lg = 2 pi r D3505 !!3506 !! Areas:3507 !! s = pi r^23508 !! s_A = pi r_A^23509 !! s_I = pi r_I^23510 !!3511 !! Life expectancy:3512 !! tau_I = 3 C* ((C*/C*t)^3/2 - 1) / r_I3513 !!3514 !! Time deratives:3515 !! dD/dt = B - (D-A)/tau_I - 2 Lg C* D3516 !! dA/dt = B - A/tau_A + 2 Lg_I C* (D-A) - 2 Lg_A C* A3517 !! 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 a03519 !!3520 3521 DO i = 1, klon3522 IF (wk_adv(i)) THEN3523 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)**23532 as_wk(i) = 3.14*arad_wk(i)**23533 is_wk(i) = 3.14*irad_wk(i)**23534 !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 ENDIF3539 ENDDO3540 3541 3542 DO i = 1, klon3543 IF (wk_adv(i)) THEN3544 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_cv3547 3548 d_sig_gen(i) = wgen(i)*aa03549 d_sig_death(i) = - isigmaw(i)*tau_wk_inv_min3550 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)*dtimesub3554 d_sig_death(i) = d_sig_death(i)*dtimesub3555 d_sig_col(i) = d_sig_col(i)*dtimesub3556 d_sig_spread(i) = d_sig_spread(i)*dtimesub3557 d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)3558 #ifdef IOPHYS_WK3559 IF (phys_sub) call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw)3560 #endif3561 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_targ3568 !! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))3569 #ifdef IOPHYS_WK3570 IF (phys_sub) THEN3571 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 ENDIF3579 #endif3580 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)*aa03583 d_asig_spread(i) = agfl(i)*cstar(i)3584 !3585 d_asig_death(i) = d_asig_death(i)*dtimesub3586 d_asig_aicol(i) = d_asig_aicol(i)*dtimesub3587 d_asig_iicol(i) = d_asig_iicol(i)*dtimesub3588 d_asig_spread(i) = d_asig_spread(i)*dtimesub3589 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_WK3591 IF (phys_sub) call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw)3592 #endif3593 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_targ3598 #ifdef IOPHYS_WK3599 IF (phys_sub) THEN3600 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 ENDIF3607 #endif3608 d_dens_gen(i) = wgen(i)3609 d_dens_death(i) = - iwdens(i)*tau_wk_inv_min3610 d_dens_col(i) = - 2.*gfl(i)*cstar(i)*wdens(i)3611 !3612 d_dens_gen(i) = d_dens_gen(i)*dtimesub3613 d_dens_death(i) = d_dens_death(i)*dtimesub3614 d_dens_col(i) = d_dens_col(i)*dtimesub3615 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_targ3621 #ifdef IOPHYS_WK3622 IF (phys_sub) THEN3623 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 ENDIF3628 #endif3629 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)*dtimesub3635 d_adens_icol(i) = d_adens_icol(i)*dtimesub3636 d_adens_acol(i) = d_adens_acol(i)*dtimesub3637 d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)3638 #ifdef IOPHYS_WK3639 IF (phys_sub) THEN3640 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 ENDIF3645 #endif3646 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_targ3650 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)*dtimesub3656 !! d_arad(i) = d_arad(i)*dtimesub3657 !! 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 ENDIF3661 ENDDO3662 3663 IF (prt_level >= 10) THEN3664 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 ENDIF3671 sigmaw=sigmaw+d_sigmaw3672 asigmaw=asigmaw+d_asigmaw3673 wdens=wdens+d_wdens3674 awdens=awdens+d_awdens3675 3676 RETURN3677 END SUBROUTINE wake_popdyn_33678 3679 2631 END MODULE lmdz_wake -
LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
r5255 r5256 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_s_a_wk,d_dens_a_wk, d_dens_wk, & ! due to wakes331 d_s_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 ISO1417 #ifdef ISOVERIF1418 write(*,*) 'physiq 1421b: qx(1,1,:)=',qx(1,1,:)1419 #endif1420 #endif1421 1416 1422 1417 … … 4450 4445 dt_a, dq_a, cv_gen, & 4451 4446 sigd, cin, & 4452 wake_deltat, wake_deltaq, wake_s, awake_ s, wake_dens, awake_dens, &4447 wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens, & 4453 4448 wake_dth, wake_h, & 4454 4449 !! wake_pe, wake_fip, wake_gfl, & … … 4460 4455 wake_omg, wake_dp_deltomg, & 4461 4456 wake_spread, wake_Cstar, d_deltat_wk_gw, & 4462 d_deltat_wk, d_deltaq_wk, d_s_wk, d_ s_a_wk, d_dens_wk, d_dens_a_wk &4457 d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk & 4463 4458 #ifdef ISO 4464 4459 & ,xt_seri,dxt_dwn,dxt_a & … … 4500 4495 4501 4496 CALL add_wake_tend & 4502 (d_deltat_wk, d_deltaq_wk, d_s_wk, d_ s_a_wk, d_dens_wk, d_dens_a_wk, wake_k, &4497 (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, & 4503 4498 'wake', abortphy & 4504 4499 #ifdef ISO … … 4775 4770 ,zqla,ztva & 4776 4771 #ifdef ISO 4777 & ,xt_ therm,d_xt_ajs &4772 & ,xt_env,d_xt_ajs & 4778 4773 #ifdef DIAGISO 4779 4774 & ,q_the,xt_the & … … 6562 6557 IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',& 6563 6558 minval(d_q_emiss),maxval(d_q_emiss) 6564 #ifdef ISO6565 abort_message='isotopes pas encore dans emission volc H2O'6566 CALL abort_physic(modname,abort_message,1)6567 #endif6568 6559 6569 6560 CALL add_phys_tend(du0, dv0, dt0, d_q_emiss, dql0, dqi0, dqbs0, paprs, & … … 6957 6948 6958 6949 print*,'avt add phys rep',abortphy 6959 #ifdef ISO 6960 abort_message='isotopes pas encore dans rep' 6961 CALL abort_physic(modname,abort_message,1) 6962 #endif 6950 6963 6951 CALL add_phys_tend & 6964 6952 (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,dqbs0,paprs,&
Note: See TracChangeset
for help on using the changeset viewer.
