MODULE lmdz_wake USE lmdz_wake_ini, ONLY: CPPKEY_IOPHYS_WK IMPLICIT NONE; PRIVATE PUBLIC wake CONTAINS SUBROUTINE wake(klon, klev, znatsurf, p, ph, pi, dtime, & tb0, qb0, omgb, & dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, & sigd_con, Cin, & deltatw, deltaqw, sigmaw, asigmaw, wdens, awdens, & ! state variables dth, hw, wape, fip, gfl, & dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, & dtke, dqke, omg, dp_deltomg, wkspread, cstar, & d_deltat_gw, & ! tendencies d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2) ! tendencies ! ************************************************************** ! * ! WAKE * ! retour a un Pupper fixe * ! * ! written by : GRANDPEIX Jean-Yves 09/03/2000 * ! modified by : ROEHRIG Romain 01/29/2007 * ! ************************************************************** USE lmdz_wake_ini, ONLY: wake_ini USE lmdz_wake_ini, ONLY: prt_level, epsim1, RG, RD USE lmdz_wake_ini, ONLY: stark, wdens_ref, coefgw, alpk, wk_pupper USE lmdz_wake_ini, ONLY: crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl USE lmdz_wake_ini, ONLY: ok_bug_gfl USE lmdz_wake_ini, ONLY: iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensinit, wdensthreshold USE lmdz_wake_ini, ONLY: sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc USE lmdz_wake_ini, ONLY: iflag_wk_profile USE lmdz_wake_ini, ONLY: smallestreal, wk_nsub IMPLICIT NONE ! ============================================================================ ! But : Decrire le comportement des poches froides apparaissant dans les ! grands systemes convectifs, et fournir l'energie disponible pour ! le declenchement de nouvelles colonnes convectives. ! State variables : ! deltatw : temperature difference between wake and off-wake regions ! deltaqw : specific humidity difference between wake and off-wake regions ! sigmaw : fractional area covered by wakes. ! asigmaw : fractional area covered by active wakes. ! wdens : number of wakes per unit area ! awdens : number of active wakes per unit area ! Variable de sortie : ! wape : WAke Potential Energy ! fip : Front Incident Power (W/m2) - ALP ! gfl : Gust Front Length per unit area (m-1) ! dtls : large scale temperature tendency due to wake ! dqls : large scale humidity tendency due to wake ! hw : wake top hight (given by hw*deltatw(1)/2=wape) ! dp_omgb : vertical gradient of large scale omega ! awdens : densite de poches actives ! wdens : densite de poches ! omgbdth: flux of Delta_Theta transported by LS omega ! dtKE : differential heating (wake - unpertubed) ! dqKE : differential moistening (wake - unpertubed) ! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) ! dp_deltomg : vertical gradient of omg (s-1) ! wkspread : spreading term in d_t_wake and d_q_wake ! deltatw : updated temperature difference (T_w-T_u). ! deltaqw : updated humidity difference (q_w-q_u). ! sigmaw : updated wake fractional area. ! asigmaw : updated active wake fractional area. ! d_deltat_gw : delta T tendency due to GW ! Variables d'entree : ! aire : aire de la maille ! tb0 : horizontal average of temperature (K) ! qb0 : horizontal average of humidity (kg/kg) ! omgb : vitesse verticale moyenne sur la maille (Pa/s) ! dtdwn: source de chaleur due aux descentes (K/s) ! dqdwn: source d'humidite due aux descentes (kg/kg/s) ! dta : source de chaleur due courants satures et detrain (K/s) ! dqa : source d'humidite due aux courants satures et detra (kg/kg/s) ! wgen : number of wakes generated per unit area and per sec (/m^2/s) ! amdwn: flux de masse total des descentes, par unite de ! surface de la maille (kg/m2/s) ! amup : flux de masse total des ascendances, par unite de ! surface de la maille (kg/m2/s) ! sigd_con: ! Cin : convective inhibition ! p : pressions aux milieux des couches (Pa) ! ph : pressions aux interfaces (Pa) ! pi : (p/p_0)**kapa (adim) ! dtime: increment temporel (s) ! Variables internes : ! rho : mean density at P levels ! rhoh : mean density at Ph levels ! tb : mean temperature | may change within ! qb : mean humidity | sub-time-stepping ! thb : mean potential temperature ! thx : potential temperature in (x) area ! tx : temperature in (x) area ! qx : humidity in (x) area ! dp_omgb: vertical gradient og LS omega ! omgbw : wake average vertical omega ! dp_omgbw: vertical gradient of omgbw ! omgbdq : flux of Delta_q transported by LS omega ! dth : potential temperature diff. wake-undist. ! th1 : first pot. temp. for vertical advection (=thx) ! th2 : second pot. temp. for vertical advection (=thw) ! q1 : first humidity for vertical advection ! q2 : second humidity for vertical advection ! d_deltatw : redistribution term for deltatw ! d_deltaqw : redistribution term for deltaqw ! deltatw0 : initial deltatw ! deltaqw0 : initial deltaqw ! hw0 : wake top hight (defined as the altitude at which deltatw=0) ! amflux : horizontal mass flux through wake boundary ! wdens_ref: initial number of wakes per unit area (3D) or per ! unit length (2D), at the beginning of each time step ! Tgw : 1 sur la periode de onde de gravite ! Cgw : vitesse de propagation de onde de gravite ! LL : distance between 2 wakes ! ------------------------------------------------------------------------- ! Declaration de variables ! ------------------------------------------------------------------------- ! Arguments en entree ! -------------------- INTEGER, INTENT(IN) :: klon, klev INTEGER, DIMENSION (klon), INTENT(IN) :: znatsurf REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi REAL, DIMENSION (klon, klev + 1), INTENT(IN) :: ph REAL, DIMENSION (klon, klev), INTENT(IN) :: omgb REAL, INTENT(IN) :: dtime REAL, DIMENSION (klon, klev), INTENT(IN) :: tb0, qb0 REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa REAL, DIMENSION (klon), INTENT(IN) :: wgen REAL, DIMENSION (klon), INTENT(IN) :: sigd_con REAL, DIMENSION (klon), INTENT(IN) :: Cin ! Input/Output ! State variables REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw REAL, DIMENSION (klon), INTENT(INOUT) :: asigmaw REAL, DIMENSION (klon), INTENT(INOUT) :: wdens REAL, DIMENSION (klon), INTENT(INOUT) :: awdens ! Sorties ! -------- REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth REAL, DIMENSION (klon, klev), INTENT(OUT) :: tx, qx REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke REAL, DIMENSION (klon, klev), INTENT(OUT) :: wkspread ! unused (jyg) REAL, DIMENSION (klon, klev), INTENT(OUT) :: omgbdth, omg REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields ! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens) REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2 ! Variables internes ! ------------------- ! Variables a fixer REAL :: delta_t_min REAL :: dtimesub REAL :: wdens0 ! IM 080208 LOGICAL, DIMENSION (klon) :: gwake ! Variables de sauvegarde REAL, DIMENSION (klon, klev) :: deltatw0 REAL, DIMENSION (klon, klev) :: deltaqw0 REAL, DIMENSION (klon, klev) :: tb, qb ! Variables liees a la dynamique de population 1 REAL, DIMENSION(klon) :: act REAL, DIMENSION(klon) :: rad_wk, tau_wk_inv REAL, DIMENSION(klon) :: f_shear REAL, DIMENSION(klon) :: drdt ! Variables liees a la dynamique de population 2 REAL, DIMENSION(klon) :: cont_fact ! Variables liees a la dynamique de population 3 REAL, DIMENSION(klon) :: arad_wk, irad_wk !! REAL, DIMENSION(klon) :: d_sig_gen, d_sig_death, d_sig_col REAL, DIMENSION(klon) :: wape1_act, wape2_act LOGICAL, DIMENSION (klon) :: kill_wake REAL :: drdt_pos REAL :: tau_wk_inv_min ! Some components of the tendencies of state variables REAL, DIMENSION (klon) :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2 REAL, DIMENSION (klon) :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2 REAL, DIMENSION (klon) :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2 REAL, DIMENSION (klon) :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2 ! Variables pour les GW REAL, DIMENSION (klon) :: ll REAL, DIMENSION (klon, klev) :: n2 REAL, DIMENSION (klon, klev) :: cgw REAL, DIMENSION (klon, klev) :: tgw ! Variables liees au calcul de hw REAL, DIMENSION (klon) :: ptop REAL, DIMENSION (klon) :: sum_dth REAL, DIMENSION (klon) :: dthmin REAL, DIMENSION (klon) :: z, dz, hw0 INTEGER, DIMENSION (klon) :: ktop, kupper ! Variables liees au test de la forme triangulaire du profil de Delta_theta REAL, DIMENSION (klon) :: sum_half_dth REAL, DIMENSION (klon) :: dz_half ! Sub-timestep tendencies and related variables REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw REAL, DIMENSION (klon, klev) :: d_tb, d_qb REAL, DIMENSION (klon) :: d_wdens, d_awdens, d_sigmaw, d_asigmaw REAL, DIMENSION (klon) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd REAL, DIMENSION (klon) :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd REAL, DIMENSION (klon) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd REAL, DIMENSION (klon) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd REAL, DIMENSION (klon) :: agfl !! gust front length of active wakes !! per unit area REAL, DIMENSION (klon) :: alpha, alpha_tot REAL, DIMENSION (klon) :: q0_min, q1_min LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw ! Autres variables internes INTEGER :: isubstep, k, i, igout REAL :: wdensmin REAL :: sigmaw_targ REAL :: wdens_targ REAL :: d_sigmaw_targ REAL :: d_wdens_targ REAL, DIMENSION (klon) :: sum_thx, sum_tx, sum_qx, sum_thvx REAL, DIMENSION (klon) :: sum_dq REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn REAL, DIMENSION (klon) :: av_thx, av_tx, av_qx, av_thvx REAL, DIMENSION (klon) :: av_dth, av_dq REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn REAL, DIMENSION (klon, klev) :: rho REAL, DIMENSION (klon, klev + 1) :: rhoh REAL, DIMENSION (klon, klev) :: zh REAL, DIMENSION (klon, klev + 1) :: zhh REAL, DIMENSION (klon, klev) :: thb, thx REAL, DIMENSION (klon, klev) :: omgbw REAL, DIMENSION (klon) :: pupper REAL, DIMENSION (klon) :: omgtop REAL, DIMENSION (klon, klev) :: dp_omgbw REAL, DIMENSION (klon) :: ztop, dztop REAL, DIMENSION (klon, klev) :: alpha_up REAL, DIMENSION (klon) :: rre1, rre2 REAL :: rrd1, rrd2 REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2 REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq REAL, DIMENSION (klon, klev) :: omgbdq REAL, DIMENSION (klon) :: ff, gg REAL, DIMENSION (klon) :: wape2, cstar2, heff REAL, DIMENSION (klon, klev) :: crep REAL, DIMENSION (klon, klev) :: ppi ! cc nrlmd REAL, DIMENSION (klon) :: death_rate !! REAL, DIMENSION (klon) :: nat_rate REAL, DIMENSION (klon, klev) :: entr REAL, DIMENSION (klon, klev) :: detr REAL, DIMENSION(klon) :: sigmaw_in, asigmaw_in ! pour les prints REAL, DIMENSION(klon) :: wdens_in, awdens_in ! pour les prints !!! LOGICAL :: phys_sub=.TRUE. LOGICAL :: phys_sub = .FALSE. LOGICAL :: first_call = .TRUE. !!-- variables liees au nouveau calcul de ptop et hw REAL, DIMENSION (klon, klev) :: int_dth REAL, DIMENSION (klon, klev) :: zzz, dzzz REAL :: epsil REAL, DIMENSION (klon) :: ptop1 INTEGER, DIMENSION (klon) :: ktop1 REAL, DIMENSION (klon) :: omega REAL, DIMENSION (klon) :: h_zzz !PRINT*,'WAKE LJYFz' ! ------------------------------------------------------------------------- ! Initialisations ! ------------------------------------------------------------------------- ! ALON = 3.e5 ! alon = 1.E6 ! Provisionnal; to be suppressed when f_shear is parameterized f_shear(:) = 1. ! 0. for strong shear, 1. for weak shear ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) ! coefgw : Coefficient pour les ondes de gravite ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) ! wdens : Densite surfacique de poche froide ! ------------------------------------------------------------------------- ! cc nrlmd coefgw=10 ! coefgw=1 ! wdens0 = 1.0/(alon**2) ! cc nrlmd wdens = 1.0/(alon**2) ! cc nrlmd stark = 0.50 ! CRtest ! cc nrlmd alpk=0.1 ! alpk = 1.0 ! alpk = 0.5 ! alpk = 0.05 igout = klon / 2 + 1 / klon ! sub-time-stepping parameters dtimesub = dtime / wk_nsub IF (first_call) THEN IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) THEN CALL iophys_ini(dtimesub) ELSE CALL iophys_ini(dtime) ENDIF END IF first_call = .FALSE. ENDIF !(first_call) IF (iflag_wk_pop_dyn == 0) THEN ! Initialisation de toutes des densites a wdens_ref. ! Les densites peuvent evoluer si les poches debordent ! (voir au tout debut de la boucle sur les substeps) !jyg< !! wdens(:) = wdens_ref DO i = 1, klon wdens(i) = wdens_ref(znatsurf(i) + 1) ENDDO !>jyg ENDIF ! (iflag_wk_pop_dyn == 0) IF (iflag_wk_pop_dyn >=1) THEN IF (iflag_wk_pop_dyn == 3) THEN wdensmin = wdensthreshold ELSE wdensmin = wdensinit ENDIF ENDIF ! PRINT*,'stark',stark ! PRINT*,'alpk',alpk ! PRINT*,'wdens',wdens ! PRINT*,'coefgw',coefgw ! cc ! Minimum value for |T_wake - T_undist|. Used for wake top definition ! ------------------------------------------------------------------------- delta_t_min = 0.2 ! 1. - Save initial values, initialize tendencies, initialize output fields ! ------------------------------------------------------------------------ !jyg< !! DO k = 1, klev !! DO i = 1, klon !! ppi(i, k) = pi(i, k) !! deltatw0(i, k) = deltatw(i, k) !! deltaqw0(i, k) = deltaqw(i, k) !! tb(i, k) = tb0(i, k) !! qb(i, k) = qb0(i, k) !! dtls(i, k) = 0. !! dqls(i, k) = 0. !! d_deltat_gw(i, k) = 0. !! d_tb(i, k) = 0. !! d_qb(i, k) = 0. !! d_deltatw(i, k) = 0. !! d_deltaqw(i, k) = 0. !! ! IM 060508 beg !! d_deltatw2(i, k) = 0. !! d_deltaqw2(i, k) = 0. !! ! IM 060508 end !! END DO !! END DO ppi(:, :) = pi(:, :) deltatw0(:, :) = deltatw(:, :) deltaqw0(:, :) = deltaqw(:, :) tb(:, :) = tb0(:, :) qb(:, :) = qb0(:, :) dtls(:, :) = 0. dqls(:, :) = 0. d_deltat_gw(:, :) = 0. d_tb(:, :) = 0. d_qb(:, :) = 0. d_deltatw(:, :) = 0. d_deltaqw(:, :) = 0. d_deltatw2(:, :) = 0. d_deltaqw2(:, :) = 0. d_sig_gen2(:) = 0. d_sig_death2(:) = 0. d_sig_col2(:) = 0. d_sig_spread2(:) = 0. d_asig_death2(:) = 0. d_asig_iicol2(:) = 0. d_asig_aicol2(:) = 0. d_asig_spread2(:) = 0. d_asig_bnd2(:) = 0. d_asigmaw2(:) = 0. d_dens_gen2(:) = 0. d_dens_death2(:) = 0. d_dens_col2(:) = 0. d_dens_bnd2(:) = 0. d_wdens2(:) = 0. d_adens_bnd2(:) = 0. d_awdens2(:) = 0. d_adens_death2(:) = 0. d_adens_icol2(:) = 0. d_adens_acol2(:) = 0. IF (iflag_wk_act == 0) THEN act(:) = 0. ELSEIF (iflag_wk_act == 1) THEN act(:) = 1. ENDIF !! DO i = 1, klon !! sigmaw_in(i) = sigmaw(i) !! END DO sigmaw_in(:) = sigmaw(:) asigmaw_in(:) = asigmaw(:) !>jyg IF (iflag_wk_pop_dyn >= 1) THEN awdens_in(:) = awdens(:) wdens_in(:) = wdens(:) !! wdens(:) = wdens(:) + wgen(:)*dtime !! d_wdens2(:) = wgen(:)*dtime !! ELSE ENDIF ! (iflag_wk_pop_dyn >= 1) ! sigmaw1=sigmaw ! IF (sigd_con.GT.sigmaw1) THEN ! PRINT*, 'sigmaw,sigd_con', sigmaw, sigd_con ! ENDIF IF (iflag_wk_pop_dyn >= 1) THEN DO i = 1, klon d_dens_gen2(i) = 0. d_dens_death2(i) = 0. d_dens_col2(i) = 0. d_awdens2(i) = 0. IF (wdens(i) < wdensthreshold) THEN !! wdens_targ = max(wdens(i),wdensmin) wdens_targ = max(wdens(i), wdensinit) d_dens_bnd2(i) = wdens_targ - wdens(i) d_wdens2(i) = wdens_targ - wdens(i) wdens(i) = wdens_targ ELSE d_dens_bnd2(i) = 0. d_wdens2(i) = 0. ENDIF !! (wdens(i) < wdensthreshold) END DO IF (iflag_wk_pop_dyn >= 2) THEN DO i = 1, klon IF (awdens(i) < wdensthreshold) THEN !! wdens_targ = min(max(awdens(i),wdensmin),wdens(i)) wdens_targ = min(max(awdens(i), wdensinit), wdens(i)) d_adens_bnd2(i) = wdens_targ - awdens(i) d_awdens2(i) = wdens_targ - awdens(i) awdens(i) = wdens_targ ELSE wdens_targ = min(awdens(i), wdens(i)) d_adens_bnd2(i) = wdens_targ - awdens(i) d_awdens2(i) = wdens_targ - awdens(i) awdens(i) = wdens_targ ENDIF END DO ENDIF ! (iflag_wk_pop_dyn >= 2) ELSE DO i = 1, klon d_awdens2(i) = 0. d_wdens2(i) = 0. END DO ENDIF ! (iflag_wk_pop_dyn >= 1) DO i = 1, klon sigmaw_targ = min(max(sigmaw(i), sigmad), 0.99) d_sig_bnd2(i) = sigmaw_targ - sigmaw(i) d_sigmaw2(i) = sigmaw_targ - sigmaw(i) sigmaw(i) = sigmaw_targ END DO IF (iflag_wk_pop_dyn == 3) THEN DO i = 1, klon IF ((wdens(i) - awdens(i)) <= smallestreal) THEN sigmaw_targ = sigmaw(i) ELSE sigmaw_targ = min(max(asigmaw(i), sigmad), sigmaw(i)) ENDIF d_asig_bnd2(i) = sigmaw_targ - asigmaw(i) d_asigmaw2(i) = sigmaw_targ - asigmaw(i) asigmaw(i) = sigmaw_targ END DO ENDIF ! (iflag_wk_pop_dyn == 3) wape(:) = 0. wape2(:) = 0. d_sigmaw(:) = 0. d_asigmaw(:) = 0. ktopw(:) = 0 !jyg IF (prt_level>=10) THEN PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout) PRINT *, 'wake-1, deltatw(igout,k) ', (k, deltatw(igout, k), k = 1, klev) PRINT *, 'wake-1, deltaqw(igout,k) ', (k, deltaqw(igout, k), k = 1, klev) PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k, amdwn(igout, k), k = 1, klev) PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k, dtdwn(igout, k), k = 1, klev) PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k, dqdwn(igout, k), k = 1, klev) PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k, amup(igout, k), k = 1, klev) PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k, dta(igout, k), k = 1, klev) PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k, dqa(igout, k), k = 1, klev) ENDIF ! 2. - Prognostic part ! -------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon z(i) = 0. ktop(i) = 0 kupper(i) = 0 sum_thx(i) = 0. sum_tx(i) = 0. sum_qx(i) = 0. sum_thvx(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thx(i) = 0. av_tx(i) = 0. av_qx(i) = 0. av_thvx(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END DO ! Distance between wakes DO i = 1, klon ll(i) = (1 - sqrt(sigmaw(i))) / sqrt(wdens(i)) END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! WRITE(*,*)'wake 1',i,k,RD,tb(i,k) rho(i, k) = p(i, k) / (RD * tb(i, k)) ! WRITE(*,*)'wake 2',rho(i,k) IF (k==1) THEN ! WRITE(*,*)'wake 3',i,k,rd,tb(i,k) rhoh(i, k) = ph(i, k) / (RD * tb(i, k)) ! WRITE(*,*)'wake 4',i,k,rd,tb(i,k) zhh(i, k) = 0 ELSE ! WRITE(*,*)'wake 5',rd,(tb(i,k)+tb(i,k-1)) rhoh(i, k) = ph(i, k) * 2. / (RD * (tb(i, k) + tb(i, k - 1))) ! WRITE(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1) zhh(i, k) = (ph(i, k) - ph(i, k - 1)) / (-rhoh(i, k) * RG) + zhh(i, k - 1) END IF ! WRITE(*,*)'wake 7',ppi(i,k) thb(i, k) = tb(i, k) / ppi(i, k) thx(i, k) = (tb(i, k) - deltatw(i, k) * sigmaw(i)) / ppi(i, k) tx(i, k) = tb(i, k) - deltatw(i, k) * sigmaw(i) qx(i, k) = qb(i, k) - deltaqw(i, k) * sigmaw(i) ! WRITE(*,*)'wake 8',(RD*(tb(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k) / ppi(i, k) END DO END DO DO k = 1, klev - 1 DO i = 1, klon IF (k==1) THEN n2(i, k) = 0 ELSE n2(i, k) = amax1(0., -RG**2 / thb(i, k) * rho(i, k) * (thb(i, k + 1) - thb(i, k - 1)) / & (p(i, k + 1) - p(i, k - 1))) END IF zh(i, k) = (zhh(i, k) + zhh(i, k + 1)) / 2 cgw(i, k) = sqrt(n2(i, k)) * zh(i, k) tgw(i, k) = coefgw * cgw(i, k) / ll(i) END DO END DO DO i = 1, klon n2(i, klev) = 0 zh(i, klev) = 0 cgw(i, klev) = 0 tgw(i, klev) = 0 END DO ! Choose an integration bound well above wake top ! ----------------------------------------------------------------- ! Determine Wake top pressure (Ptop) from buoyancy integral ! -------------------------------------------------------- Do i = 1, klon wk_adv(i) = .True. Enddo Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, & dth, hw0, rho, delta_t_min, & ktop, wk_adv, h_zzz, ptop1, ktop1) !!print'("pkupper APPEL ",7i6)',0,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper IF (prt_level>=10) THEN PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) ENDIF ! -5/ Set deltatw & deltaqw to 0 above kupper DO k = 1, klev DO i = 1, klon IF (k>=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. d_deltatw2(i, k) = -deltatw0(i, k) d_deltaqw2(i, k) = -deltaqw0(i, k) END IF END DO END DO ! Vertical gradient of LS omega DO k = 1, klev DO i = 1, klon IF (k<=kupper(i)) THEN dp_omgb(i, k) = (omgb(i, k + 1) - omgb(i, k)) / (ph(i, k + 1) - ph(i, k)) END IF END DO END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvx to 1st level virt. pot. temp. DO i = 1, klon z(i) = 1. dz(i) = 1. sum_thvx(i) = thx(i, 1) * (1. + epsim1 * qx(i, 1)) * dz(i) sum_dth(i) = 0. END DO DO k = 1, klev DO i = 1, klon dz(i) = -(amax1(ph(i, k + 1), ptop(i)) - ph(i, k)) / (rho(i, k) * RG) IF (dz(i)>0) THEN ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau z(i) = z(i) + dz(i) sum_thx(i) = sum_thx(i) + thx(i, k) * dz(i) sum_tx(i) = sum_tx(i) + tx(i, k) * dz(i) sum_qx(i) = sum_qx(i) + qx(i, k) * dz(i) sum_thvx(i) = sum_thvx(i) + thx(i, k) * (1. + epsim1 * qx(i, k)) * dz(i) sum_dth(i) = sum_dth(i) + dth(i, k) * dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k) * dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k) * dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k) * dz(i) END IF END DO END DO DO i = 1, klon hw0(i) = z(i) END DO ! 2.1 - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon av_thx(i) = sum_thx(i) / hw0(i) av_tx(i) = sum_tx(i) / hw0(i) av_qx(i) = sum_qx(i) / hw0(i) av_thvx(i) = sum_thvx(i) / hw0(i) ! av_thve = sum_thve/hw0 av_dth(i) = sum_dth(i) / hw0(i) av_dq(i) = sum_dq(i) / hw0(i) av_dtdwn(i) = sum_dtdwn(i) / hw0(i) av_dqdwn(i) = sum_dqdwn(i) / hw0(i) wape(i) = -RG * hw0(i) * (av_dth(i) + & epsim1 * (av_thx(i) * av_dq(i) + av_dth(i) * av_qx(i) + av_dth(i) * av_dq(i))) / av_thvx(i) END DO IF (CPPKEY_IOPHYS_WK) THEN IF (.NOT.phys_sub) CALL iophys_ecrit('wape_a', 1, 'wape_a', 'J/kg', wape) END IF ! 2.2 Prognostic variable update ! ------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. d_deltatw2(i, k) = -deltatw0(i, k) d_deltaqw2(i, k) = -deltaqw0(i, k) END IF END DO END DO DO i = 1, klon IF (wape(i)<0.) THEN !! sigmaw(i) = amax1(sigmad, sigd_con(i)) sigmaw_targ = max(sigmad, sigd_con(i)) d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) sigmaw(i) = sigmaw_targ ENDIF !! (wape(i)<0.) ENDDO IF (iflag_wk_pop_dyn == 3) THEN DO i = 1, klon IF (wape(i)<0.) THEN sigmaw_targ = max(sigmad, sigd_con(i)) d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) asigmaw(i) = sigmaw_targ ENDIF !! (wape(i)<0.) ENDDO ENDIF !! (iflag_wk_pop_dyn == 3) DO i = 1, klon IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin fip(i) = 0. gwake(i) = .FALSE. ELSE hw(i) = hw0(i) cstar(i) = stark * sqrt(2. * wape(i)) gwake(i) = .TRUE. END IF END DO ! Check qx and qw positivity ! -------------------------- DO i = 1, klon q0_min(i) = min((qb(i, 1) - sigmaw(i) * deltaqw(i, 1)), & (qb(i, 1) + (1. - sigmaw(i)) * deltaqw(i, 1))) END DO DO k = 2, klev DO i = 1, klon q1_min(i) = min((qb(i, k) - sigmaw(i) * deltaqw(i, k)), & (qb(i, k) + (1. - sigmaw(i)) * deltaqw(i, k))) IF (q1_min(i)<=q0_min(i)) THEN q0_min(i) = q1_min(i) END IF END DO END DO DO i = 1, klon ok_qx_qw(i) = q0_min(i) >= 0. alpha(i) = 1. alpha_tot(i) = 1. END DO IF (prt_level>=10) THEN PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', & sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ENDIF ! C ----------------------------------------------------------------- ! Sub-time-stepping ! ----------------- ! wk_nsub and dtimesub definitions moved to begining of routine. !! wk_nsub = 10 !! dtimesub = dtime/wk_nsub ! ------------------------------------------------------------------------ ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ------------------------------------------------------------------------ DO isubstep = 1, wk_nsub ! ------------------------------------------------------------------------ ! wk_adv is the LOGICAL flag enabling wake evolution in the time advance ! loop DO i = 1, klon wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1. END DO IF (prt_level>=10) THEN PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', & isubstep, wk_adv(igout), cstar(igout), wape(igout), ptop(igout) ENDIF ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement ! negatif de ktop a kupper -------- ! cc On calcule pour cela une densite wdens0 pour laquelle on ! aurait un entrainement nul --- !jyg< ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou ! les poches sont insuffisantes pour accueillir tout le flux de masse ! des descentes unsaturees. Nous faisons alors l'hypothese que la ! convection profonde cree directement de nouvelles poches, sans passer ! par les thermiques. La nouvelle valeur de wdens est alors imposee. DO i = 1, klon ! c PRINT *,' isubstep,wk_adv(i),cstar(i),wape(i) ', ! c $ isubstep,wk_adv(i),cstar(i),wape(i) IF (wk_adv(i) .AND. cstar(i)>0.01) THEN IF (iflag_wk_profile == 0) THEN omg(i, kupper(i) + 1) = -RG * amdwn(i, kupper(i) + 1) / sigmaw(i) + & RG * amup(i, kupper(i) + 1) / (1. - sigmaw(i)) ELSE omg(i, kupper(i) + 1) = 0. ENDIF wdens0 = (sigmaw(i) / (4. * 3.14)) * & ((1. - sigmaw(i)) * omg(i, kupper(i) + 1) / ((ph(i, 1) - pupper(i)) * cstar(i)))**(2) IF (prt_level >= 10) THEN PRINT*, 'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', & omg(i, kupper(i) + 1), wdens0, wdens(i), cstar(i), ph(i, 1) - pupper(i) ENDIF IF (wdens(i)<=wdens0 * 1.1) THEN IF (iflag_wk_pop_dyn >= 1) THEN d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i) d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i) ENDIF wdens(i) = wdens0 END IF END IF END DO IF (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) THEN !!-------------------------------------------------------- !!Bug : computing gfl and rad_wk before changing sigmaw !! This bug exists only for iflag_wk_pop_dyn=0. Otherwise, gfl and rad_wk !! are computed within wake_popdyn !!-------------------------------------------------------- DO i = 1, klon IF (wk_adv(i)) THEN gfl(i) = 2. * sqrt(3.14 * wdens(i) * sigmaw(i)) rad_wk(i) = sqrt(sigmaw(i) / (3.14 * wdens(i))) END IF END DO ENDIF ! (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) !!-------------------------------------------------------- DO i = 1, klon IF (wk_adv(i)) THEN sigmaw_targ = min(sigmaw(i), sigmaw_max) d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) sigmaw(i) = sigmaw_targ END IF END DO IF (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) THEN !!-------------------------------------------------------- !!Fix : computing gfl and rad_wk after changing sigmaw !!-------------------------------------------------------- DO i = 1, klon IF (wk_adv(i)) THEN gfl(i) = 2. * sqrt(3.14 * wdens(i) * sigmaw(i)) rad_wk(i) = sqrt(sigmaw(i) / (3.14 * wdens(i))) END IF END DO ENDIF ! (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) !!-------------------------------------------------------- IF (iflag_wk_pop_dyn >= 1) THEN ! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0. ! Here, it has to be set to zero. death_rate(:) = 0. ENDIF IF (iflag_wk_pop_dyn >= 3) THEN DO i = 1, klon IF (wk_adv(i)) THEN sigmaw_targ = min(asigmaw(i), sigmaw_max) d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) asigmaw(i) = sigmaw_targ ENDIF ENDDO ENDIF !!-------------------------------------------------------- !!-------------------------------------------------------- IF (iflag_wk_pop_dyn == 1) THEN CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, & wdensmin, & dtimesub, gfl, rad_wk, f_shear, drdt_pos, & d_awdens, d_wdens, d_sigmaw, & iflag_wk_act, wk_adv, cin, wape, & drdt, & d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & d_wdens_targ, d_sigmaw_targ) !!-------------------------------------------------------- ELSEIF (iflag_wk_pop_dyn == 2) THEN CALL wake_popdyn_2 (klon, klev, wk_adv, dtimesub, wgen, & wdensmin, & sigmaw, wdens, awdens, & !! state variables gfl, cstar, cin, wape, rad_wk, & d_sigmaw, d_wdens, d_awdens, & !! tendencies cont_fact, & d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd) sigmaw = sigmaw - d_sigmaw wdens = wdens - d_wdens awdens = awdens - d_awdens !!-------------------------------------------------------- ELSEIF (iflag_wk_pop_dyn == 3) THEN IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) THEN CALL iophys_ecrit('ptop', 1, 'ptop', 'Pa', ptop) CALL iophys_ecrit('sigmaw', 1, 'sigmaw', '', sigmaw) CALL iophys_ecrit('asigmaw', 1, 'asigmaw', '', asigmaw) CALL iophys_ecrit('wdens', 1, 'wdens', '1/m2', wdens) CALL iophys_ecrit('awdens', 1, 'awdens', '1/m2', awdens) CALL iophys_ecrit('rad_wk', 1, 'rad_wk', 'm', rad_wk) CALL iophys_ecrit('arad_wk', 1, 'arad_wk', 'm', arad_wk) CALL iophys_ecrit('irad_wk', 1, 'irad_wk', 'm', irad_wk) ENDIF END IF CALL wake_popdyn_3 (klon, klev, phys_sub, wk_adv, dtimesub, wgen, & wdensmin, & sigmaw, asigmaw, wdens, awdens, & !! state variables gfl, agfl, cstar, cin, wape, & rad_wk, arad_wk, irad_wk, & d_sigmaw, d_asigmaw, d_wdens, d_awdens, & !! tendencies d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, & d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd) sigmaw = sigmaw - d_sigmaw asigmaw = asigmaw - d_asigmaw wdens = wdens - d_wdens awdens = awdens - d_awdens !!-------------------------------------------------------- ELSEIF (iflag_wk_pop_dyn == 0) THEN ! cc nrlmd DO i = 1, klon IF (wk_adv(i)) THEN ! cc nrlmd Introduction du taux de mortalite des poches et ! test sur sigmaw_max=0.4 ! cc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub IF (sigmaw(i)>=sigmaw_max) THEN death_rate(i) = gfl(i) * cstar(i) / sigmaw(i) ELSE death_rate(i) = 0. END IF d_sigmaw(i) = gfl(i) * cstar(i) * dtimesub - death_rate(i) * sigmaw(i) * & dtimesub ! $ - nat_rate(i)*sigmaw(i)*dtimesub ! c PRINT*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i)', ! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i) ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub ! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!! ! wdens = wdens0/(10.*sigmaw) ! sigmaw =max(sigmaw,sigd_con) ! sigmaw =max(sigmaw,sigmad) END IF END DO ENDIF ! (iflag_wk_pop_dyn == 1) ... ELSEIF (iflag_wk_pop_dyn == 0) !!-------------------------------------------------------- !!-------------------------------------------------------- IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) THEN CALL iophys_ecrit('wdensa', 1, 'wdensa', 'm', wdens) CALL iophys_ecrit('awdensa', 1, 'awdensa', 'm', awdens) CALL iophys_ecrit('sigmawa', 1, 'sigmawa', 'm', sigmaw) CALL iophys_ecrit('asigmawa', 1, 'asigmawa', 'm', asigmaw) ENDIF END IF ! calcul de la difference de vitesse verticale poche - zone non perturbee ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit ! IM 060208 au niveau k=1... !JYG 161013 Correction : maintenant omg est dimensionne a klev. DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dp_deltomg(i, k) = 0. END IF END DO END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd omg(i, k) = 0. END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN z(i) = 0. omg(i, 1) = 0. dp_deltomg(i, 1) = -(gfl(i) * cstar(i)) / (sigmaw(i) * (1 - sigmaw(i))) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN dz(i) = -(ph(i, k) - ph(i, k - 1)) / (rho(i, k - 1) * RG) z(i) = z(i) + dz(i) dp_deltomg(i, k) = dp_deltomg(i, 1) omg(i, k) = dp_deltomg(i, 1) * z(i) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN dztop(i) = -(ptop(i) - ph(i, ktop(i))) / (rho(i, ktop(i)) * RG) ztop(i) = z(i) + dztop(i) omgtop(i) = dp_deltomg(i, 1) * ztop(i) END IF END DO IF (prt_level>=10) THEN PRINT *, 'wake-4.2, omg(igout,k) ', (k, omg(igout, k), k = 1, klev) PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', & omgtop(igout), ptop(igout), ktop(igout) ENDIF ! ----------------- ! From m/s to Pa/s ! ----------------- DO i = 1, klon IF (wk_adv(i)) THEN omgtop(i) = -rho(i, ktop(i)) * RG * omgtop(i) !! LJYF dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) dp_deltomg(i, 1) = omgtop(i) / min(ptop(i) - ph(i, 1), -smallestreal) END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN omg(i, k) = -rho(i, k) * RG * omg(i, k) dp_deltomg(i, k) = dp_deltomg(i, 1) END IF END DO END DO ! raccordement lineaire de omg de ptop a pupper DO i = 1, klon IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN IF (iflag_wk_profile == 0) THEN omg(i, kupper(i) + 1) = -RG * amdwn(i, kupper(i) + 1) / sigmaw(i) + & RG * amup(i, kupper(i) + 1) / (1. - sigmaw(i)) ELSE omg(i, kupper(i) + 1) = 0. ENDIF dp_deltomg(i, kupper(i)) = (omgtop(i) - omg(i, kupper(i) + 1)) / & (ptop(i) - pupper(i)) END IF END DO ! c DO i=1,klon ! c PRINT*,'Pente entre 0 et kupper (reference)' ! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) ! c PRINT*,'Pente entre ktop et kupper' ! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i)) ! c ENDDO ! c DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN dp_deltomg(i, k) = dp_deltomg(i, kupper(i)) omg(i, k) = omgtop(i) + (ph(i, k) - ptop(i)) * dp_deltomg(i, kupper(i)) END IF END DO END DO !! PRINT *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev) ! cc nrlmd ! c DO i=1,klon ! c PRINT*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1) ! c END DO ! cc ! -- Compute wake average vertical velocity omgbw DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN omgbw(i, k) = omgb(i, k) + (1. - sigmaw(i)) * omg(i, k) END IF END DO END DO ! -- and its vertical gradient dp_omgbw DO k = 1, klev - 1 DO i = 1, klon IF (wk_adv(i)) THEN dp_omgbw(i, k) = (omgbw(i, k + 1) - omgbw(i, k)) / (ph(i, k + 1) - ph(i, k)) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN dp_omgbw(i, klev) = 0. END IF END DO ! -- Upstream coefficients for omgb velocity ! -- (alpha_up(k) is the coefficient of the value at level k) ! -- (1-alpha_up(k) is the coefficient of the value at level k-1) DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN alpha_up(i, k) = 0. IF (omgb(i, k)>0.) alpha_up(i, k) = 1. END IF END DO END DO ! Matrix expressing [The,deltatw] from [Th1,Th2] DO i = 1, klon IF (wk_adv(i)) THEN rre1(i) = 1. - sigmaw(i) rre2(i) = sigmaw(i) END IF END DO rrd1 = -1. rrd2 = 1. ! -- Get [Th1,Th2], dth and [q1,q2] DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i) + 1) THEN dth(i, k) = deltatw(i, k) / ppi(i, k) th1(i, k) = thb(i, k) - sigmaw(i) * dth(i, k) ! undisturbed area th2(i, k) = thb(i, k) + (1. - sigmaw(i)) * dth(i, k) ! wake q1(i, k) = qb(i, k) - sigmaw(i) * deltaqw(i, k) ! undisturbed area q2(i, k) = qb(i, k) + (1. - sigmaw(i)) * deltaqw(i, k) ! wake END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd d_th1(i, 1) = 0. d_th2(i, 1) = 0. d_dth(i, 1) = 0. d_q1(i, 1) = 0. d_q2(i, 1) = 0. d_dq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i) + 1) THEN d_th1(i, k) = th1(i, k - 1) - th1(i, k) d_th2(i, k) = th2(i, k - 1) - th2(i, k) d_dth(i, k) = dth(i, k - 1) - dth(i, k) d_q1(i, k) = q1(i, k - 1) - q1(i, k) d_q2(i, k) = q2(i, k - 1) - q2(i, k) d_dq(i, k) = deltaqw(i, k - 1) - deltaqw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN omgbdth(i, 1) = 0. omgbdq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i) + 1) THEN ! loop on interfaces omgbdth(i, k) = omgb(i, k) * (dth(i, k - 1) - dth(i, k)) omgbdq(i, k) = omgb(i, k) * (deltaqw(i, k - 1) - deltaqw(i, k)) END IF END DO END DO !! IF (prt_level>=10) THEN IF (prt_level>=10 .AND. wk_adv(igout)) THEN PRINT *, 'wake-4.3, th1(igout,k) ', (k, th1(igout, k), k = 1, kupper(igout)) PRINT *, 'wake-4.3, th2(igout,k) ', (k, th2(igout, k), k = 1, kupper(igout)) PRINT *, 'wake-4.3, dth(igout,k) ', (k, dth(igout, k), k = 1, kupper(igout)) PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k, omgbdth(igout, k), k = 1, kupper(igout)) ENDIF ! ----------------------------------------------------------------- DO k = 1, klev - 1 DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i) - 1) THEN ! ----------------------------------------------------------------- ! Compute redistribution (advective) term d_deltatw(i, k) = dtimesub / (ph(i, k) - ph(i, k + 1)) * & (rrd1 * omg(i, k) * sigmaw(i) * d_th1(i, k) - & rrd2 * omg(i, k + 1) * (1. - sigmaw(i)) * d_th2(i, k + 1) - & (1. - alpha_up(i, k)) * omgbdth(i, k) - & alpha_up(i, k + 1) * omgbdth(i, k + 1)) * ppi(i, k) ! PRINT*,'d_d,k_ptop_provis(i)eltatw=', k, d_deltatw(i,k) d_deltaqw(i, k) = dtimesub / (ph(i, k) - ph(i, k + 1)) * & (rrd1 * omg(i, k) * sigmaw(i) * d_q1(i, k) - & rrd2 * omg(i, k + 1) * (1. - sigmaw(i)) * d_q2(i, k + 1) - & (1. - alpha_up(i, k)) * omgbdq(i, k) - & alpha_up(i, k + 1) * omgbdq(i, k + 1)) ! PRINT*,'d_deltaqw=', k, d_deltaqw(i,k) ! and increment large scale tendencies ! C ! ----------------------------------------------------------------- d_tb(i, k) = dtimesub * ((rre1(i) * omg(i, k) * sigmaw(i) * d_th1(i, k) - & rre2(i) * omg(i, k + 1) * (1. - sigmaw(i)) * d_th2(i, k + 1)) / & (ph(i, k) - ph(i, k + 1)) & - sigmaw(i) * (1. - sigmaw(i)) * dth(i, k) * (omg(i, k) - omg(i, k + 1)) / & (ph(i, k) - ph(i, k + 1))) * ppi(i, k) d_qb(i, k) = dtimesub * ((rre1(i) * omg(i, k) * sigmaw(i) * d_q1(i, k) - & rre2(i) * omg(i, k + 1) * (1. - sigmaw(i)) * d_q2(i, k + 1)) / & (ph(i, k) - ph(i, k + 1)) & - sigmaw(i) * (1. - sigmaw(i)) * deltaqw(i, k) * (omg(i, k) - omg(i, k + 1)) / & (ph(i, k) - ph(i, k + 1))) ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN d_tb(i, k) = dtimesub * (rre1(i) * omg(i, k) * sigmaw(i) * d_th1(i, k) / (ph(i, k) - ph(i, k + 1))) * ppi(i, k) d_qb(i, k) = dtimesub * (rre1(i) * omg(i, k) * sigmaw(i) * d_q1(i, k) / (ph(i, k) - ph(i, k + 1))) END IF ! cc END DO END DO ! ------------------------------------------------------------------ IF (prt_level>=10) THEN PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k, d_deltatw(igout, k), k = 1, klev) PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k, d_deltaqw(igout, k), k = 1, klev) ENDIF ! Increment state variables !jyg< IF (iflag_wk_pop_dyn >= 1) THEN DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN detr(i, k) = - d_sig_death(i) - d_sig_col(i) entr(i, k) = d_sig_gen(i) ENDIF ENDDO ENDDO ELSE ! (iflag_wk_pop_dyn >= 1) DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN detr(i, k) = 0. entr(i, k) = 0. ENDIF ENDDO ENDDO ENDIF ! (iflag_wk_pop_dyn >= 1) DO k = 1, klev DO i = 1, klon ! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN IF (wk_adv(i) .AND. k<=kupper(i)) THEN ! cc ! Coefficient de repartition crep(i, k) = crep_sol * (ph(i, kupper(i)) - ph(i, k)) / & (ph(i, kupper(i)) - ph(i, 1)) crep(i, k) = crep(i, k) + crep_upper * (ph(i, 1) - ph(i, k)) / & (ph(i, 1) - ph(i, kupper(i))) ! Reintroduce compensating subsidence term. ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) ! . /(1-sigmaw) ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) ! . /(1-sigmaw) dtke(i, k) = (dtdwn(i, k) / sigmaw(i) - dta(i, k) / (1. - sigmaw(i))) dqke(i, k) = (dqdwn(i, k) / sigmaw(i) - dqa(i, k) / (1. - sigmaw(i))) ! PRINT*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k) ! cc nrlmd Prise en compte du taux de mortalite ! cc Definitions de entr, detr !jyg< !! detr(i, k) = 0. !! !! entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + & !! sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) !! entr(i, k) = entr(i, k) + gfl(i) * cstar(i) + & sigmaw(i) * (1. - sigmaw(i)) * dp_deltomg(i, k) !>jyg wkspread(i, k) = (entr(i, k) - detr(i, k)) / sigmaw(i) ! cc wkspread(i,k) = ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ ! cc $ sigmaw(i) ! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU ! Jingmei ! WRITE(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k), ! & Tgw(i,k),deltatw(i,k) d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k) * deltatw(i, k) * & dtimesub ! WRITE(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k) ff(i) = d_deltatw(i, k) / dtimesub ! Sans GW ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-wkspread(k)*deltatw(k)) ! GW formule 1 ! deltatw(k) = deltatw(k)+dtimesub* ! $ (ff+dtKE(k) - wkspread(k)*deltatw(k)-Tgw(k)*deltatw(k)) ! GW formule 2 IF (dtimesub * tgw(i, k)<1.E-10) THEN d_deltatw(i, k) = dtimesub * (ff(i) + dtke(i, k) - & entr(i, k) * deltatw(i, k) / sigmaw(i) - & (death_rate(i) * sigmaw(i) + detr(i, k)) * deltatw(i, k) / (1. - sigmaw(i)) - & ! cc tgw(i, k) * deltatw(i, k)) ELSE d_deltatw(i, k) = 1 / tgw(i, k) * (1 - exp(-dtimesub * tgw(i, k))) * & (ff(i) + dtke(i, k) - & entr(i, k) * deltatw(i, k) / sigmaw(i) - & (death_rate(i) * sigmaw(i) + detr(i, k)) * deltatw(i, k) / (1. - sigmaw(i)) - & tgw(i, k) * deltatw(i, k)) END IF dth(i, k) = deltatw(i, k) / ppi(i, k) gg(i) = d_deltaqw(i, k) / dtimesub d_deltaqw(i, k) = dtimesub * (gg(i) + dqke(i, k) - & entr(i, k) * deltaqw(i, k) / sigmaw(i) - & (death_rate(i) * sigmaw(i) + detr(i, k)) * deltaqw(i, k) / (1. - sigmaw(i))) ! cc ! cc nrlmd ! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) ! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) ! cc END IF END DO END DO ! Scale tendencies so that water vapour remains positive in w and x. CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qb, d_qb, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! Alpha_tot = Product of all the alpha's DO i = 1, klon IF (wk_adv(i)) THEN alpha_tot(i) = alpha_tot(i) * alpha(i) END IF END DO ! cc nrlmd ! c PRINT*,'alpha' ! c do i=1,klon ! c PRINT*,alpha(i) ! c END DO ! cc DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN d_tb(i, k) = alpha(i) * d_tb(i, k) d_qb(i, k) = alpha(i) * d_qb(i, k) d_deltatw(i, k) = alpha(i) * d_deltatw(i, k) d_deltaqw(i, k) = alpha(i) * d_deltaqw(i, k) d_deltat_gw(i, k) = alpha(i) * d_deltat_gw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN d_sigmaw(i) = alpha(i) * d_sigmaw(i) END IF END DO ! Update large scale variables and wake variables ! IM 060208 manque DO i + remplace DO k=1,kupper(i) ! IM 060208 DO k = 1,kupper(i) DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN dtls(i, k) = dtls(i, k) + d_tb(i, k) dqls(i, k) = dqls(i, k) + d_qb(i, k) ! cc nrlmd d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k) d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k) ! cc END IF END DO END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN tb(i, k) = tb0(i, k) + dtls(i, k) qb(i, k) = qb0(i, k) + dqls(i, k) thb(i, k) = tb(i, k) / ppi(i, k) deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k) deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k) dth(i, k) = deltatw(i, k) / ppi(i, k) ! c PRINT*,'k,qx,qw',k,qb(i,k)-sigmaw(i)*deltaqw(i,k) ! c $ ,qb(i,k)+(1-sigmaw(i))*deltaqw(i,k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN sigmaw(i) = sigmaw(i) + d_sigmaw(i) d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i) END IF END DO !jyg< IF (iflag_wk_pop_dyn >= 1) THEN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sigmaw !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Cumulatives DO i = 1, klon IF (wk_adv(i)) THEN d_sig_gen2(i) = d_sig_gen2(i) + d_sig_gen(i) d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i) d_sig_col2(i) = d_sig_col2(i) + d_sig_col(i) d_sig_spread2(i) = d_sig_spread2(i) + d_sig_spread(i) d_sig_bnd2(i) = d_sig_bnd2(i) + d_sig_bnd(i) END IF END DO ! Bounds DO i = 1, klon IF (wk_adv(i)) THEN sigmaw_targ = max(sigmaw(i), sigmad) d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) sigmaw(i) = sigmaw_targ END IF END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Cumulatives DO i = 1, klon IF (wk_adv(i)) THEN wdens(i) = wdens(i) + d_wdens(i) d_wdens2(i) = d_wdens2(i) + d_wdens(i) d_dens_gen2(i) = d_dens_gen2(i) + d_dens_gen(i) d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i) d_dens_col2(i) = d_dens_col2(i) + d_dens_col(i) d_dens_bnd2(i) = d_dens_bnd2(i) + d_dens_bnd(i) END IF END DO ! Bounds DO i = 1, klon IF (wk_adv(i)) THEN wdens_targ = max(wdens(i), wdensmin) d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i) d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i) wdens(i) = wdens_targ END IF END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Cumulatives DO i = 1, klon IF (wk_adv(i)) THEN awdens(i) = awdens(i) + d_awdens(i) d_awdens2(i) = d_awdens2(i) + d_awdens(i) END IF END DO ! Bounds DO i = 1, klon IF (wk_adv(i)) THEN wdens_targ = min(max(awdens(i), 0.), wdens(i)) d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i) awdens(i) = wdens_targ END IF END DO IF (iflag_wk_pop_dyn >= 2) THEN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn >= 2!!!!!! ! Cumulatives DO i = 1, klon IF (wk_adv(i)) THEN d_adens_death2(i) = d_adens_death2(i) + d_adens_death(i) d_adens_icol2(i) = d_adens_icol2(i) + d_adens_icol(i) d_adens_acol2(i) = d_adens_acol2(i) + d_adens_acol(i) d_adens_bnd2(i) = d_adens_bnd2(i) + d_adens_bnd(i) END IF END DO ! Bounds DO i = 1, klon IF (wk_adv(i)) THEN wdens_targ = min(max(awdens(i), 0.), wdens(i)) d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) awdens(i) = wdens_targ END IF END DO IF (iflag_wk_pop_dyn == 3) THEN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! asigmaw for iflag_wk_pop_dyn = 3!!!!!! ! Cumulatives DO i = 1, klon IF (wk_adv(i)) THEN asigmaw(i) = asigmaw(i) + d_asigmaw(i) d_asigmaw2(i) = d_asigmaw2(i) + d_asigmaw(i) d_asig_death2(i) = d_asig_death2(i) + d_asig_death(i) d_asig_spread2(i) = d_asig_spread2(i) + d_asig_spread(i) d_asig_iicol2(i) = d_asig_iicol2(i) + d_asig_iicol(i) d_asig_aicol2(i) = d_asig_aicol2(i) + d_asig_aicol(i) d_asig_bnd2(i) = d_asig_bnd2(i) + d_asig_bnd(i) END IF END DO ! Bounds DO i = 1, klon IF (wk_adv(i)) THEN ! asigmaw lower bound set to sigmad/2 in order to allow asigmaw values lower than sigmad. !! sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i)) sigmaw_targ = min(max(asigmaw(i), sigmad / 2.), sigmaw(i)) d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) asigmaw(i) = sigmaw_targ END IF END DO IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) THEN CALL iophys_ecrit('wdensb', 1, 'wdensb', 'm', wdens) CALL iophys_ecrit('awdensb', 1, 'awdensb', 'm', awdens) CALL iophys_ecrit('sigmawb', 1, 'sigmawb', 'm', sigmaw) CALL iophys_ecrit('asigmawb', 1, 'asigmawb', 'm', asigmaw) CALL iophys_ecrit('d_wdens2', 1, 'd_wdens2', '', d_wdens2) CALL iophys_ecrit('d_dens_gen2', 1, 'd_dens_gen2', '', d_dens_gen2) CALL iophys_ecrit('d_dens_death2', 1, 'd_dens_death2', '', d_dens_death2) CALL iophys_ecrit('d_dens_col2', 1, 'd_dens_col2', '', d_dens_col2) CALL iophys_ecrit('d_dens_bnd2', 1, 'd_dens_bnd2', '', d_dens_bnd2) CALL iophys_ecrit('d_awdens2', 1, 'd_awdens2', '', d_awdens2) CALL iophys_ecrit('d_adens_death2', 1, 'd_adens_death2', '', d_adens_death2) CALL iophys_ecrit('d_adens_icol2', 1, 'd_adens_icol2', '', d_adens_icol2) CALL iophys_ecrit('d_adens_acol2', 1, 'd_adens_acol2', '', d_adens_acol2) CALL iophys_ecrit('d_adens_bnd2', 1, 'd_adens_bnd2', '', d_adens_bnd2) CALL iophys_ecrit('d_sigmaw2', 1, 'd_sigmaw2', '', d_sigmaw2) CALL iophys_ecrit('d_sig_gen2', 1, 'd_sig_gen2', 'm', d_sig_gen2) CALL iophys_ecrit('d_sig_spread2', 1, 'd_sig_spread2', '', d_sig_spread2) CALL iophys_ecrit('d_sig_col2', 1, 'd_sig_col2', '', d_sig_col2) CALL iophys_ecrit('d_sig_death2', 1, 'd_sig_death2', '', d_sig_death2) CALL iophys_ecrit('d_sig_bnd2', 1, 'd_sig_bnd2', '', d_sig_bnd2) CALL iophys_ecrit('d_asigmaw2', 1, 'd_asigmaw2', '', d_asigmaw2) CALL iophys_ecrit('d_asig_spread2', 1, 'd_asig_spread2', 'm', d_asig_spread2) CALL iophys_ecrit('d_asig_aicol2', 1, 'd_asig_aicol2', 'm', d_asig_aicol2) CALL iophys_ecrit('d_asig_iicol2', 1, 'd_asig_iicol2', 'm', d_asig_iicol2) CALL iophys_ecrit('d_asig_death2', 1, 'd_asig_death2', 'm', d_asig_death2) CALL iophys_ecrit('d_asig_bnd2', 1, 'd_asig_bnd2', 'm', d_asig_bnd2) ENDIF END IF ENDIF ! (iflag_wk_pop_dyn == 3) ENDIF ! (iflag_wk_pop_dyn >= 2) ENDIF ! (iflag_wk_pop_dyn >= 1) Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, & dth, hw, rho, delta_t_min, & ktop, wk_adv, h_zzz, ptop1, ktop1) !! print'("pkupper APPEL ",7i6)',isubstep,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper ! 5/ Set deltatw & deltaqw to 0 above kupper DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. d_deltatw2(i, k) = -deltatw0(i, k) d_deltaqw2(i, k) = -deltaqw0(i, k) END IF END DO END DO ! -------------Cstar computation--------------------------------- DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_thx(i) = 0. sum_tx(i) = 0. sum_qx(i) = 0. sum_thvx(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thx(i) = 0. av_tx(i) = 0. av_qx(i) = 0. av_thvx(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvx to 1st level virt. pot. temp. DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd z(i) = 1. dz(i) = 1. sum_thvx(i) = thx(i, 1) * (1. + epsim1 * qx(i, 1)) * dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dz(i) = -(max(ph(i, k + 1), ptop(i)) - ph(i, k)) / (rho(i, k) * RG) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thx(i) = sum_thx(i) + thx(i, k) * dz(i) sum_tx(i) = sum_tx(i) + tx(i, k) * dz(i) sum_qx(i) = sum_qx(i) + qx(i, k) * dz(i) sum_thvx(i) = sum_thvx(i) + thx(i, k) * (1. + epsim1 * qx(i, k)) * dz(i) sum_dth(i) = sum_dth(i) + dth(i, k) * dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k) * dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k) * dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k) * dz(i) END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd av_thx(i) = sum_thx(i) / hw0(i) av_tx(i) = sum_tx(i) / hw0(i) av_qx(i) = sum_qx(i) / hw0(i) av_thvx(i) = sum_thvx(i) / hw0(i) av_dth(i) = sum_dth(i) / hw0(i) av_dq(i) = sum_dq(i) / hw0(i) av_dtdwn(i) = sum_dtdwn(i) / hw0(i) av_dqdwn(i) = sum_dqdwn(i) / hw0(i) wape(i) = -RG * hw0(i) * (av_dth(i) + epsim1 * (av_thx(i) * av_dq(i) + & av_dth(i) * av_qx(i) + av_dth(i) * av_dq(i))) / av_thvx(i) END IF END DO ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. d_deltatw2(i, k) = -deltatw0(i, k) d_deltaqw2(i, k) = -deltaqw0(i, k) END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin !jyg< !! sigmaw(i) = max(sigmad, sigd_con(i)) sigmaw_targ = max(sigmad, sigd_con(i)) d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) sigmaw(i) = sigmaw_targ d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) asigmaw(i) = sigmaw_targ !>jyg fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark * sqrt(2. * wape(i)) gwake(i) = .TRUE. END IF END IF END DO ! ------------------------------------------------------------------------ END DO ! isubstep end sub-timestep loop ! ------------------------------------------------------------------------ ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ------------------------------------------------------------------------ IF (CPPKEY_IOPHYS_WK) THEN IF (.NOT.phys_sub) CALL iophys_ecrit('wape_b', 1, 'wape_b', 'J/kg', wape) END IF IF (prt_level>=10) THEN PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', & sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ENDIF ! ---------------------------------------------------------- ! Determine wake final state; recompute wape, cstar, ktop; ! filter out bad wakes. ! ---------------------------------------------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon ! cc nrlmd if (wk_adv(i)) then !!! nrlmd IF (ok_qx_qw(i)) THEN ! cc z(i) = 0. sum_thx(i) = 0. sum_tx(i) = 0. sum_qx(i) = 0. sum_thvx(i) = 0. sum_dth(i) = 0. sum_half_dth(i) = 0. sum_dq(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thx(i) = 0. av_tx(i) = 0. av_qx(i) = 0. av_thvx(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. dthmin(i) = -delta_t_min END IF END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc rho(i, k) = p(i, k) / (RD * tb(i, k)) IF (k==1) THEN rhoh(i, k) = ph(i, k) / (RD * tb(i, k)) zhh(i, k) = 0 ELSE rhoh(i, k) = ph(i, k) * 2. / (RD * (tb(i, k) + tb(i, k - 1))) zhh(i, k) = (ph(i, k) - ph(i, k - 1)) / (-rhoh(i, k) * RG) + zhh(i, k - 1) END IF thb(i, k) = tb(i, k) / ppi(i, k) thx(i, k) = (tb(i, k) - deltatw(i, k) * sigmaw(i)) / ppi(i, k) tx(i, k) = tb(i, k) - deltatw(i, k) * sigmaw(i) qx(i, k) = qb(i, k) - deltaqw(i, k) * sigmaw(i) dth(i, k) = deltatw(i, k) / ppi(i, k) END IF END DO END DO ! Integrals (and wake top level number) ! ----------------------------------------------------------- ! Initialize sum_thvx to 1st level virt. pot. temp. DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc z(i) = 1. dz(i) = 1. dz_half(i) = 1. sum_thvx(i) = thx(i, 1) * (1. + epsim1 * qx(i, 1)) * dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc dz(i) = -(amax1(ph(i, k + 1), ptop(i)) - ph(i, k)) / (rho(i, k) * RG) dz_half(i) = -(amax1(ph(i, k + 1), 0.5 * (ptop(i) + ph(i, 1))) - ph(i, k)) / (rho(i, k) * RG) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thx(i) = sum_thx(i) + thx(i, k) * dz(i) sum_tx(i) = sum_tx(i) + tx(i, k) * dz(i) sum_qx(i) = sum_qx(i) + qx(i, k) * dz(i) sum_thvx(i) = sum_thvx(i) + thx(i, k) * (1. + epsim1 * qx(i, k)) * dz(i) sum_dth(i) = sum_dth(i) + dth(i, k) * dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k) * dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k) * dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k) * dz(i) dthmin(i) = min(dthmin(i), dth(i, k)) END IF IF (dz_half(i)>0) THEN sum_half_dth(i) = sum_half_dth(i) + dth(i, k) * dz_half(i) END IF END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! ------------------------------------------------------------- ! Means DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc av_thx(i) = sum_thx(i) / hw0(i) av_tx(i) = sum_tx(i) / hw0(i) av_qx(i) = sum_qx(i) / hw0(i) av_thvx(i) = sum_thvx(i) / hw0(i) av_dth(i) = sum_dth(i) / hw0(i) av_dq(i) = sum_dq(i) / hw0(i) av_dtdwn(i) = sum_dtdwn(i) / hw0(i) av_dqdwn(i) = sum_dqdwn(i) / hw0(i) wape2(i) = -RG * hw0(i) * (av_dth(i) + epsim1 * (av_thx(i) * av_dq(i) + & av_dth(i) * av_qx(i) + av_dth(i) * av_dq(i))) / av_thvx(i) END IF END DO IF (CPPKEY_IOPHYS_WK) THEN IF (.NOT.phys_sub) CALL iophys_ecrit('wape2_a', 1, 'wape2_a', 'J/kg', wape2) END IF ! Prognostic variable update ! ------------------------------------------------------------ ! Filter out bad wakes IF (iflag_wk_check_trgl>=1) THEN ! Check triangular shape of dth profile DO i = 1, klon IF (ok_qx_qw(i)) THEN !! PRINT *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i) !! PRINT *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', & !! 2.*sum_dth(i)/(hw0(i)*dthmin(i)) !! PRINT *,'wake, sum_half_dth(i), sum_dth(i) ', & !! sum_half_dth(i), sum_dth(i) IF ((hw0(i) < 1.) .OR. (dthmin(i) >= -delta_t_min)) THEN wape2(i) = -1. !! PRINT *,'wake, rej 1' ELSE IF (iflag_wk_check_trgl==1.AND.abs(2. * sum_dth(i) / (hw0(i) * dthmin(i)) - 1.) > 0.5) THEN wape2(i) = -1. !! PRINT *,'wake, rej 2' ELSE IF (abs(sum_half_dth(i)) < 0.5 * abs(sum_dth(i))) THEN wape2(i) = -1. !! PRINT *,'wake, rej 3' END IF END IF END DO END IF IF (CPPKEY_IOPHYS_WK) THEN IF (.NOT.phys_sub) CALL iophys_ecrit('wape2_b', 1, 'wape2_b', 'J/kg', wape2) END IF DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN ! cc deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. d_deltatw2(i, k) = -deltatw0(i, k) d_deltaqw2(i, k) = -deltaqw0(i, k) END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (wape2(i)<0.) THEN wape2(i) = 0. cstar2(i) = 0. hw(i) = hwmin !jyg< !! sigmaw(i) = amax1(sigmad, sigd_con(i)) sigmaw_targ = max(sigmad, sigd_con(i)) d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) sigmaw(i) = sigmaw_targ d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i) asigmaw(i) = sigmaw_targ !>jyg fip(i) = 0. gwake(i) = .FALSE. ELSE IF (prt_level>=10) PRINT *, 'wape2>0' cstar2(i) = stark * sqrt(2. * wape2(i)) gwake(i) = .TRUE. END IF IF (CPPKEY_IOPHYS_WK) THEN IF (.NOT.phys_sub) CALL iophys_ecrit('cstar2', 1, 'cstar2', 'J/kg', cstar2) END IF END IF ! (ok_qx_qw(i)) END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc ktopw(i) = ktop(i) END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (ktopw(i)>0 .AND. gwake(i)) THEN ! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) ! cc heff = 600. ! Utilisation de la hauteur hw ! c heff = 0.7*hw heff(i) = hw(i) fip(i) = 0.5 * rho(i, ktopw(i)) * cstar2(i)**3 * heff(i) * 2 * & sqrt(sigmaw(i) * wdens(i) * 3.14) fip(i) = alpk * fip(i) ! jyg2 ELSE fip(i) = 0. END IF END IF END DO IF (iflag_wk_pop_dyn >= 3) THEN IF (CPPKEY_IOPHYS_WK) THEN IF (.NOT.phys_sub) THEN CALL iophys_ecrit('fip', 1, 'fip', 'J/kg', fip) CALL iophys_ecrit('hw', 1, 'hw', 'J/kg', hw) CALL iophys_ecrit('ptop', 1, 'ptop', 'J/kg', ptop) CALL iophys_ecrit('wdens', 1, 'wdens', 'J/kg', wdens) CALL iophys_ecrit('awdens', 1, 'awdens', 'm', awdens) CALL iophys_ecrit('sigmaw', 1, 'sigmaw', 'm', sigmaw) CALL iophys_ecrit('asigmaw', 1, 'asigmaw', 'm', asigmaw) CALL iophys_ecrit('rad_wk', 1, 'rad_wk', 'J/kg', rad_wk) CALL iophys_ecrit('arad_wk', 1, 'arad_wk', 'J/kg', arad_wk) CALL iophys_ecrit('irad_wk', 1, 'irad_wk', 'J/kg', irad_wk) CALL iophys_ecrit('d_wdens2', 1, 'd_wdens2', '', d_wdens2) CALL iophys_ecrit('d_dens_gen2', 1, 'd_dens_gen2', '', d_dens_gen2) CALL iophys_ecrit('d_dens_death2', 1, 'd_dens_death2', '', d_dens_death2) CALL iophys_ecrit('d_dens_col2', 1, 'd_dens_col2', '', d_dens_col2) CALL iophys_ecrit('d_dens_bnd2', 1, 'd_dens_bnd2', '', d_dens_bnd2) CALL iophys_ecrit('d_awdens2', 1, 'd_awdens2', '', d_awdens2) CALL iophys_ecrit('d_adens_death2', 1, 'd_adens_death2', '', d_adens_death2) CALL iophys_ecrit('d_adens_icol2', 1, 'd_adens_icol2', '', d_adens_icol2) CALL iophys_ecrit('d_adens_acol2', 1, 'd_adens_acol2', '', d_adens_acol2) CALL iophys_ecrit('d_adens_bnd2', 1, 'd_adens_bnd2', '', d_adens_bnd2) CALL iophys_ecrit('d_sigmaw2', 1, 'd_sigmaw2', '', d_sigmaw2) CALL iophys_ecrit('d_sig_gen2', 1, 'd_sig_gen2', 'm', d_sig_gen2) CALL iophys_ecrit('d_sig_spread2', 1, 'd_sig_spread2', '', d_sig_spread2) CALL iophys_ecrit('d_sig_col2', 1, 'd_sig_col2', '', d_sig_col2) CALL iophys_ecrit('d_sig_death2', 1, 'd_sig_death2', '', d_sig_death2) CALL iophys_ecrit('d_sig_bnd2', 1, 'd_sig_bnd2', '', d_sig_bnd2) CALL iophys_ecrit('d_asigmaw2', 1, 'd_asigmaw2', '', d_asigmaw2) CALL iophys_ecrit('d_asig_spread2', 1, 'd_asig_spread2', 'm', d_asig_spread2) CALL iophys_ecrit('d_asig_aicol2', 1, 'd_asig_aicol2', 'm', d_asig_aicol2) CALL iophys_ecrit('d_asig_iicol2', 1, 'd_asig_iicol2', 'm', d_asig_iicol2) CALL iophys_ecrit('d_asig_death2', 1, 'd_asig_death2', 'm', d_asig_death2) CALL iophys_ecrit('d_asig_bnd2', 1, 'd_asig_bnd2', 'm', d_asig_bnd2) ENDIF ! (.NOT.phys_sub) END IF ENDIF ! (iflag_wk_pop_dyn >= 3) ! Limitation de sigmaw ! cc nrlmd ! DO i=1,klon ! IF (OK_qx_qw(i)) THEN ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max ! ENDIF ! ENDDO ! cc !jyg< IF (iflag_wk_pop_dyn >= 1) THEN DO i = 1, klon kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold) !! .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin) ENDDO ELSE ! (iflag_wk_pop_dyn >= 1) DO i = 1, klon kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i) ENDDO ENDIF ! (iflag_wk_pop_dyn >= 1) !>jyg DO k = 1, klev DO i = 1, klon !!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & !!jyg .NOT. ok_qx_qw(i)) THEN IF (kill_wake(i)) THEN ! cc dtls(i, k) = 0. dqls(i, k) = 0. deltatw(i, k) = 0. deltaqw(i, k) = 0. d_deltatw2(i, k) = -deltatw0(i, k) d_deltaqw2(i, k) = -deltaqw0(i, k) END IF ! (kill_wake(i)) END DO END DO DO i = 1, klon !!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & !!jyg .NOT. ok_qx_qw(i)) THEN IF (kill_wake(i)) THEN ktopw(i) = 0 wape(i) = 0. cstar(i) = 0. !!jyg Outside SUBROUTINE "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes !! hw(i) = hwmin !jyg !! sigmaw(i) = sigmad !jyg hw(i) = 0. !jyg fip(i) = 0. !! sigmaw(i) = 0. !jyg sigmaw_targ = 0. d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) !! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i) ! _in = correction jyg 20220124 sigmaw(i) = sigmaw_targ IF (iflag_wk_pop_dyn >= 3) THEN sigmaw_targ = 0. d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i) !! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) d_asigmaw2(i) = sigmaw_targ - asigmaw_in(i) ! _in = correction jyg 20220124 asigmaw(i) = sigmaw_targ ELSE asigmaw(i) = 0. ENDIF ! (iflag_wk_pop_dyn >= 3) IF (iflag_wk_pop_dyn >= 1) THEN !! awdens(i) = 0. !! wdens(i) = 0. wdens_targ = 0. d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i) !! d_wdens2(i) = wdens_targ - wdens(i) d_wdens2(i) = wdens_targ - wdens_in(i) ! jyg 20220916 wdens(i) = wdens_targ wdens_targ = 0. !!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens. IF (iflag_wk_pop_dyn >= 2) THEN d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) ENDIF ! (iflag_wk_pop_dyn >= 2) !! d_awdens2(i) = wdens_targ - awdens(i) d_awdens2(i) = wdens_targ - awdens_in(i) ! jyg 20220916 awdens(i) = wdens_targ !! IF (iflag_wk_pop_dyn == 2) THEN !! d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) !! ENDIF ! (iflag_wk_pop_dyn == 2) ENDIF ! (iflag_wk_pop_dyn >= 1) ELSE ! (kill_wake(i)) wape(i) = wape2(i) cstar(i) = cstar2(i) END IF ! (kill_wake(i)) ! c PRINT*,'wape wape2 ktopw OK_qx_qw =', ! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) END DO IF (prt_level>=10) THEN PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', & wape(igout), wape2(igout), ktopw(igout), OK_qx_qw(igout) ENDIF IF (CPPKEY_IOPHYS_WK) THEN IF (.NOT.phys_sub) CALL iophys_ecrit('wape_c', 1, 'wape_c', 'J/kg', wape) END IF ! ----------------------------------------------------------------- ! Get back to tendencies per second DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN !jyg< !! IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN IF (ok_qx_qw(i)) THEN !>jyg ! cc dtls(i, k) = dtls(i, k) / dtime dqls(i, k) = dqls(i, k) / dtime d_deltatw2(i, k) = d_deltatw2(i, k) / dtime d_deltaqw2(i, k) = d_deltaqw2(i, k) / dtime d_deltat_gw(i, k) = d_deltat_gw(i, k) / dtime ! c PRINT*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k) ! c $ ,death_rate(i)*sigmaw(i) END IF END DO END DO !jyg< IF (iflag_wk_pop_dyn >= 1) THEN DO i = 1, klon IF (ok_qx_qw(i)) THEN d_sig_gen2(i) = d_sig_gen2(i) / dtime d_sig_death2(i) = d_sig_death2(i) / dtime d_sig_col2(i) = d_sig_col2(i) / dtime d_sig_spread2(i) = d_sig_spread2(i) / dtime d_sig_bnd2(i) = d_sig_bnd2(i) / dtime d_sigmaw2(i) = d_sigmaw2(i) / dtime d_dens_gen2(i) = d_dens_gen2(i) / dtime d_dens_death2(i) = d_dens_death2(i) / dtime d_dens_col2(i) = d_dens_col2(i) / dtime d_dens_bnd2(i) = d_dens_bnd2(i) / dtime d_awdens2(i) = d_awdens2(i) / dtime d_wdens2(i) = d_wdens2(i) / dtime ENDIF ENDDO IF (iflag_wk_pop_dyn >= 2) THEN DO i = 1, klon IF (ok_qx_qw(i)) THEN d_adens_death2(i) = d_adens_death2(i) / dtime d_adens_icol2(i) = d_adens_icol2(i) / dtime d_adens_acol2(i) = d_adens_acol2(i) / dtime d_adens_bnd2(i) = d_adens_bnd2(i) / dtime ENDIF ENDDO IF (iflag_wk_pop_dyn == 3) THEN DO i = 1, klon IF (ok_qx_qw(i)) THEN d_asig_death2(i) = d_asig_death2(i) / dtime d_asig_iicol2(i) = d_asig_iicol2(i) / dtime d_asig_aicol2(i) = d_asig_aicol2(i) / dtime d_asig_spread2(i) = d_asig_spread2(i) / dtime d_asig_bnd2(i) = d_asig_bnd2(i) / dtime ENDIF ENDDO ENDIF ! (iflag_wk_pop_dyn == 3) ENDIF ! (iflag_wk_pop_dyn >= 2) ENDIF ! (iflag_wk_pop_dyn >= 1) !>jyg END SUBROUTINE wake SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qb, d_qb, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! ------------------------------------------------------ ! Dtermination du coefficient alpha tel que les tendances ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent ! a une humidite positive dans la zone (x) et dans la zone (w). ! ------------------------------------------------------ IMPLICIT NONE ! Input REAL qb(nlon, nl), d_qb(nlon, nl) REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl) REAL sigmaw(nlon), d_sigmaw(nlon) LOGICAL wk_adv(nlon) INTEGER nl, nlon ! Output REAL alpha(nlon) ! Internal variables REAL zeta(nlon, nl) REAL alpha1(nlon) REAL x, a, b, c, discrim REAL epsilon_loc INTEGER i, k DO k = 1, nl DO i = 1, nlon IF (wk_adv(i)) THEN IF ((deltaqw(i, k) + d_deltaqw(i, k))>=0.) THEN zeta(i, k) = 0. ELSE zeta(i, k) = 1. END IF END IF END DO DO i = 1, nlon IF (wk_adv(i)) THEN x = qb(i, k) + (zeta(i, k) - sigmaw(i)) * deltaqw(i, k) + d_qb(i, k) + & (zeta(i, k) - sigmaw(i)) * d_deltaqw(i, k) - d_sigmaw(i) * & (deltaqw(i, k) + d_deltaqw(i, k)) a = -d_sigmaw(i) * d_deltaqw(i, k) b = d_qb(i, k) + (zeta(i, k) - sigmaw(i)) * d_deltaqw(i, k) - & deltaqw(i, k) * d_sigmaw(i) c = qb(i, k) + (zeta(i, k) - sigmaw(i)) * deltaqw(i, k) + epsilon_loc discrim = b * b - 4. * a * c ! PRINT*, 'x, a, b, c, discrim', x, a, b, c, discrim IF (a + b>=0.) THEN !! Condition suffisante pour la positivite de ovap alpha1(i) = 1. ELSE IF (x>=0.) THEN alpha1(i) = 1. ELSE IF (a>0.) THEN alpha1(i) = 0.9 * min((2. * c) / (-b + sqrt(discrim)), & (-b + sqrt(discrim)) / (2. * a)) ELSE IF (a==0.) THEN alpha1(i) = 0.9 * (-c / b) ELSE ! PRINT*,'a,b,c discrim',a,b,c discrim alpha1(i) = 0.9 * max((2. * c) / (-b + sqrt(discrim)), & (-b + sqrt(discrim)) / (2. * a)) END IF END IF END IF alpha(i) = min(alpha(i), alpha1(i)) END IF END DO END DO END SUBROUTINE wake_vec_modulation SUBROUTINE pkupper(klon, klev, ptop, ph, p, pupper, kupper, & dth, hw_, rho, delta_t_min_in, & ktop, wk_adv, h_zzz, ptop1, ktop1) USE lmdz_wake_ini, ONLY: wk_pupper USE lmdz_wake_ini, ONLY: RG USE lmdz_wake_ini, ONLY: hwmin USE lmdz_wake_ini, ONLY: iflag_wk_new_ptop, wk_delta_t_min, wk_frac_int_delta_t USE lmdz_wake_ini, ONLY: wk_int_delta_t_min IMPLICIT NONE INTEGER, INTENT(IN) :: klon, klev REAL, DIMENSION (klon, klev + 1), INTENT(IN) :: ph, p REAL, DIMENSION (klon, klev + 1), INTENT(IN) :: rho LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv REAL, DIMENSION (klon, klev + 1), INTENT(IN) :: dth REAL, INTENT(IN) :: delta_t_min_in REAL, DIMENSION (klon), INTENT(OUT) :: hw_ REAL, DIMENSION (klon), INTENT(OUT) :: ptop INTEGER, DIMENSION (klon), INTENT(OUT) :: Ktop REAL, DIMENSION (klon), INTENT(OUT) :: pupper INTEGER, DIMENSION (klon), INTENT(OUT) :: kupper REAL, DIMENSION (klon), INTENT(OUT) :: h_zzz !! REAL, DIMENSION (klon), INTENT(OUT) :: Ptop1 !! INTEGER, DIMENSION (klon), INTENT(OUT) :: ktop1 !! INTEGER :: i, k LOGICAL, DIMENSION (klon) :: wk_active REAL :: delta_t_min REAL, DIMENSION (klon) :: dthmin REAL, DIMENSION (klon) :: ptop_provis, ptop_new REAL, DIMENSION (klon) :: z, dz REAL, DIMENSION (klon) :: sum_dth INTEGER, DIMENSION (klon) :: k_ptop_provis REAL, DIMENSION (klon) :: zk_ptop_provis REAL, DIMENSION (klon) :: omega !! REAL, DIMENSION (klon, klev + 1) :: int_dth !! REAL, DIMENSION (klon, klev + 1) :: dzz !! REAL, DIMENSION (klon, klev + 1) :: zzz !! REAL, DIMENSION (klon) :: frac_int_dth !! REAL :: ddd!! INTEGER, SAVE :: ipas = 0 !INTEGER, SAVE :: compte=0 ! LJYF : a priori z, dz sum_dth sont aussi des variables internes ! Les eliminer apres verification convergence numerique !compte=compte+1 !PRINT*,'compte=',compte ! Determine Ptop from buoyancy integral ! --------------------------------------- ! - 1/ Pressure of the level where dth changes sign. !PRINT*,'WAKE LJYF' IF (iflag_wk_new_ptop==0) THEN delta_t_min = delta_t_min_in else delta_t_min = wk_delta_t_min END IF DO i = 1, klon ptop_provis(i) = ph(i, 1) k_ptop_provis(i) = 1 END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. ptop_provis(i)==ph(i, 1) .AND. & ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN dth(i, k)>-delta_t_min .AND. dth(i, k - 1)<-delta_t_min) THEN ptop_provis(i) = ((dth(i, k) + delta_t_min) * p(i, k - 1) - & (dth(i, k - 1) + delta_t_min) * p(i, k)) / (dth(i, k) - dth(i, k - 1)) k_ptop_provis(i) = k END IF END DO END DO ! - 2/ dth integral DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = -(amax1(ph(i, k + 1), ptop_provis(i)) - ph(i, k)) / (rho(i, k) * RG) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k) * dz(i) dthmin(i) = amin1(dthmin(i), dth(i, k)) END IF END IF END DO END DO ! - 3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon IF (wk_adv(i)) THEN hw_(i) = 2. * sum_dth(i) / amin1(dthmin(i), -0.5) hw_(i) = amax1(hwmin, hw_(i)) END IF END DO ! - 4/ now, get Ptop DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd ktop(i) = 0 z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = amin1(-(ph(i, k + 1) - ph(i, k)) / (rho(i, k) * RG), hw_(i) - z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k) * RG * dz(i) ktop(i) = k END IF END IF END DO END DO ! 4.5/Correct ktop and ptop DO i = 1, klon ptop_new(i) = ptop(i) END DO DO k = klev, 2, -1 DO i = 1, klon ! IM v3JYG; IF (k .GE. ktop(i) IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN dth(i, k)>=-delta_t_min .AND. dth(i, k - 1)<-delta_t_min) THEN ptop_new(i) = ((dth(i, k) + delta_t_min) * p(i, k - 1) - & (dth(i, k - 1) + delta_t_min) * p(i, k)) / (dth(i, k) - dth(i, k - 1)) END IF END DO END DO DO i = 1, klon ptop(i) = ptop_new(i) END DO DO k = klev, 1, -1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (ph(i, k + 1)=10) THEN ! PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) ! ENDIF ! ----------------------------------------------------------------------- ! nouveau calcul de hw et ptop ! ----------------------------------------------------------------------- !if (iflag_wk_new_ptop>0) THEN DO i = 1, klon ptop1(i) = ph(i, 1) ktop1(i) = 1 h_zzz(i) = 0. END DO IF (iflag_wk_new_ptop/=0) THEN int_dth(1:klon, 1:klev + 1) = 0. DO i = 1, klon IF (wk_adv(i)) THEN int_dth(i, 1) = 0. END IF END DO IF (abs(iflag_wk_new_ptop) == 1) THEN DO k = 2, klev + 1 Do i = 1, klon IF (wk_adv(i)) THEN IF (k<=k_ptop_provis(i)) THEN ddd = dth(i, k - 1) * (ph(i, k - 1) - max(ptop_provis(i), ph(i, k))) !ddd=dth(i,k-1)*(ph(i,k-1) - ph(i,k)) else ddd = 0. endif int_dth(i, k) = int_dth(i, k - 1) + ddd !ELSE ! int_dth(i,k) = 0. END IF END DO END DO else k_ptop_provis(:) = klev + 1 dthmin(:) = dth(:, 1) ! calcul de l'int??grale de dT * dP jusqu'au dernier ! niveau avec dT<0. (en s'assurant qu'on a bien un ! dT negatif plus bas) DO k = 1, klev DO i = 1, klon dthmin(i) = min(dthmin(i), dth(i, k)) ddd = dth(i, k) * (ph(i, k) - ph(i, k + 1)) IF (dthmin(i)<0.) THEN IF (k>=k_ptop_provis(i)) THEN ddd = 0. ELSE IF (dth(i, k)>=0.) THEN ddd = 0. k_ptop_provis(i) = k + 1 endif endif int_dth(i, k + 1) = int_dth(i, k) + ddd ENDDO ENDDO DO i = 1, klon IF (k_ptop_provis(i)==klev + 1 .OR. .NOT. wk_adv(i)) THEN k_ptop_provis(i) = 1 endif ENDDO endif ! (abs(iflag_wk_new_ptop) == 1 ) ! PRINT*, 'xxx, int_dth', (k,int_dth(1,k),k=1,klev) ! PRINT*, 'xxx, k_ptop_provis', k_ptop_provis(1) ! On se limite ?? des poches avec integrale dT * dp < -wk_int_delta_t_min DO i = 1, klon IF (int_dth(i, k_ptop_provis(i)) > -wk_int_delta_t_min .OR. k_ptop_provis(i)==1) THEN !if (1==0) THEN wk_active(i) = .FALSE. ptop(i) = ph(i, 1) ktop(i) = 1 hw_(i) = 0. else wk_active(i) = wk_adv(i) endif enddo DO i = 1, klon IF (wk_active(i)) THEN frac_int_dth(i) = wk_frac_int_delta_t * int_dth(i, k_ptop_provis(i)) ENDIF ENDDO DO k = 1, klev DO i = 1, klon ! PRINT*,ipas,'yyy ',k,int_dth(i,k),frac_int_dth(i) IF (wk_active(i)) THEN IF (int_dth(i, k)>=frac_int_dth(i)) THEN ktop1(i) = min(k, k_ptop_provis(i)) !ktop1(i) = k !PRINT*,ipas,'yyy ktop1= ',ktop1 ENDIF ENDIF END DO END DO !PRINT*, 'LAMINE' DO i = 1, klon IF (wk_active(i)) THEN !PRINT*, ipas,'xxx1, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ',ktop1 ddd = int_dth(i, ktop1(i) + 1) - int_dth(i, ktop1(i)) IF (ddd==0.) THEN omega(i) = 0. else omega(i) = (frac_int_dth(i) - int_dth(i, ktop1(i))) / ddd endif !! PRINT*,'OMEGA ',omega(i) END IF END DO !! PRINT*, 'xxx' DO i = 1, klon IF (wk_active(i)) THEN ! PRINT*, 'xxx, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ', & ! int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ! PRINT*, 'xxx, omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ', & !e omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ptop1(i) = min((1 - omega(i)) * ph(i, ktop1(i)) + omega(i) * ph(i, ktop1(i) + 1), ph(i, 1)) END IF END DO DO i = 1, klon IF (wk_active(i)) THEN zzz(i, 1) = 0 END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_active(i)) THEN dzz(i, k) = (ph(i, k) - ph(i, k + 1)) / (rho(i, k) * RG) zzz(i, k + 1) = zzz(i, k) + dzz(i, k) END IF END DO END DO DO i = 1, klon IF (wk_active(i)) THEN h_zzz(i) = max((1 - omega(i)) * zzz(i, ktop1(i)) + omega(i) * zzz(i, ktop1(i) + 1), hwmin) END IF END DO ENDIF ! (iflag_wk_new_ptop/=0) !if (iflag_wk_new_ptop==2) THEN IF (iflag_wk_new_ptop>0) THEN DO i = 1, klon ptop(i) = ptop1(i) ktop(i) = ktop1(i) hw_(i) = h_zzz(i) enddo !endif ENDIF kupper = 0 IF (wk_pupper<1.) THEN ! Choose an integration bound well above wake top ! ----------------------------------------------------------------- ! Pupper = 50000. ! melting level ! Pupper = 60000. ! Pupper = 80000. ! essais pour case_e DO i = 1, klon ! pupper(i) = 0.6*ph(i, 1) pupper(i) = wk_pupper * ph(i, 1) pupper(i) = max(pupper(i), 45000.) ! cc Pupper(i) = 60000. END DO ELSE DO i = 1, klon ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) ! pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-50.) pupper(i) = min(wk_pupper * ptop(i) + (1. - wk_pupper) * ph(i, 1), ptop(i) - 5000.) END DO END IF ! -5/ Determination de kupper DO k = klev, 1, -1 DO i = 1, klon IF (ph(i, k + 1)jyg+mlt !jyg< d_wdens_targ = max(d_wdens(i), wdensmin - wdens(i)) !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) d_dens_bnd(i) = d_wdens_targ - d_wdens(i) d_wdens(i) = d_wdens_targ !! d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i)) !>jyg !jyg+mlt< !! d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) & !! + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) & !! - sigmaw(i)*tau_wk_inv_min )*dtimesub d_sig_gen(i) = wgen(i) * aa0 d_sig_death(i) = - sigmaw(i) * (1. - awdens(i) / wdens(i)) * tau_wk_inv_min !! d_sig_col(i) = - 2 * f_shear(i) * sigmaw(i) * gfl(i) * drdt_pos d_sig_col(i) = - 2 * f_shear(i) * (2. * sigmaw(i) - wdens(i) * aa0) * gfl(i) * drdt_pos d_sig_spread(i) = gfl(i) * cstar(i) d_sig_gen(i) = d_sig_gen(i) * dtimesub d_sig_death(i) = d_sig_death(i) * dtimesub d_sig_col(i) = d_sig_col(i) * dtimesub d_sig_spread(i) = d_sig_spread(i) * dtimesub d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) !>jyg+mlt !jyg< d_sigmaw_targ = max(d_sigmaw(i), sigmad - sigmaw(i)) !! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) !! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) d_sigmaw(i) = d_sigmaw_targ !! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) !>jyg ENDIF ENDDO IF (prt_level >= 10) THEN PRINT *, 'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', & cstar(1), cstar(1) / cstart, rad_wk(1), tau_wk_inv(1), drdt(1) PRINT *, 'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', & wdens(1), awdens(1), act(1), d_awdens(1) PRINT *, 'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', & wgen(1), -(wdens(1) - awdens(1)) * tau_wk_inv(1), -2. * wdens(1) * gfl(1) * drdt_pos, d_wdens(1) PRINT *, 'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ENDIF END SUBROUTINE wake_popdyn_1 SUBROUTINE wake_popdyn_2(klon, klev, wk_adv, dtimesub, wgen, & wdensmin, & sigmaw, wdens, awdens, & !! states variables gfl, cstar, cin, wape, rad_wk, & d_sigmaw, d_wdens, d_awdens, & !! tendences cont_fact, & d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd) USE lmdz_wake_ini, ONLY: wake_ini USE lmdz_wake_ini, ONLY: prt_level, RG USE lmdz_wake_ini, ONLY: stark, wdens_ref USE lmdz_wake_ini, ONLY: tau_cv, rzero, aa0 !! USE lmdz_wake_ini , ONLY: iflag_wk_pop_dyn, wdensmin USE lmdz_wake_ini, ONLY: iflag_wk_pop_dyn USE lmdz_wake_ini, ONLY: sigmad, cstart, sigmaw_max IMPLICIT NONE INTEGER, INTENT(IN) :: klon, klev LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv REAL, INTENT(IN) :: dtimesub REAL, INTENT(IN) :: wdensmin REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = wake radius REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front lenght per unit area REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_wdens, d_awdens REAL, DIMENSION (klon), INTENT(OUT) :: cont_fact !! RM facteur de contact = 2 pi * rad * C* ! Some components of the tendencies of state variables REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd !! internal variables INTEGER :: i, k REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes REAL :: tau_wk_inv_min REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes REAL :: d_wdens_targ, d_sigmaw_targ !! Equations !! dD/dt = B - (D-A)/tau - f D^2 !! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^2 !! dsigma/dt = B a0 - sigma/D (D-A)/tau + Lg C* - f (D-A)^2 (sigma/D-a0) !! !! f = 2 (B (a0-sigma/D) + Lg C*) / (2 (D-A)^2 (2 sigma/D-a0) + D (1-2 sigma)) DO i = 1, klon IF (wk_adv(i)) THEN rad_wk(i) = max(sqrt(sigmaw(i) / (3.14 * wdens(i))), rzero) gfl(i) = 2. * sqrt(3.14 * wdens(i) * sigmaw(i)) END IF END DO DO i = 1, klon IF (wk_adv(i)) THEN !! tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.) tau_wk_inv(i) = max((3. * cstar(i)) / (rad_wk(i) * ((cstar(i) / cstart)**1.5 - 1)), 0.) tau_wk_inv_min = min(tau_wk_inv(i), 1. / dtimesub) tau_prime(i) = tau_cv !! cont_fact(i) = 2.*(wgen(i)*(aa0-sigmaw(i)/wdens(i)) + gfl(i)*cstar(i)) / & !! (2.*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i) - aa0) + wdens(i)*(1.-2.*sigmaw(i))) !! cont_fact(i) = 2.*3.14*rad_wk(i)*cstar(i) ! bug !! cont_fact(i) = 4.*3.14*rad_wk(i)*cstar(i) cont_fact(i) = 2. * gfl(i) * cstar(i) / wdens(i) d_sig_gen(i) = wgen(i) * aa0 d_sig_death(i) = - sigmaw(i) * (1. - awdens(i) / wdens(i)) * tau_wk_inv_min d_sig_col(i) = - cont_fact(i) * (wdens(i) - awdens(i))**2 * (2. * sigmaw(i) / wdens(i) - aa0) d_sig_spread(i) = gfl(i) * cstar(i) d_sig_gen(i) = d_sig_gen(i) * dtimesub d_sig_death(i) = d_sig_death(i) * dtimesub d_sig_col(i) = d_sig_col(i) * dtimesub d_sig_spread(i) = d_sig_spread(i) * dtimesub d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) d_sigmaw_targ = max(d_sigmaw(i), sigmad - sigmaw(i)) !! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) !! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) d_sigmaw(i) = d_sigmaw_targ !! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) d_dens_gen(i) = wgen(i) d_dens_death(i) = - (wdens(i) - awdens(i)) * tau_wk_inv_min d_dens_col(i) = - cont_fact(i) * wdens(i)**2 d_dens_gen(i) = d_dens_gen(i) * dtimesub d_dens_death(i) = d_dens_death(i) * dtimesub d_dens_col(i) = d_dens_col(i) * dtimesub d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i) d_adens_death(i) = -awdens(i) / tau_prime(i) d_adens_icol(i) = cont_fact(i) * (wdens(i) - awdens(i))**2 d_adens_acol(i) = - cont_fact(i) * awdens(i)**2 d_adens_death(i) = d_adens_death(i) * dtimesub d_adens_icol(i) = d_adens_icol(i) * dtimesub d_adens_acol(i) = d_adens_acol(i) * dtimesub d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i) !! d_wdens_targ = max(d_wdens(i), wdensmin - wdens(i)) !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) d_dens_bnd(i) = d_wdens_targ - d_wdens(i) d_wdens(i) = d_wdens_targ d_wdens_targ = min(max(d_awdens(i), -awdens(i)), wdens(i) - awdens(i)) !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) d_adens_bnd(i) = d_wdens_targ - d_awdens(i) d_awdens(i) = d_wdens_targ ENDIF ENDDO IF (prt_level >= 10) THEN PRINT *, 'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) ', & cstar(1), cstar(1) / cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) PRINT *, 'wake, wdens(1), awdens(1), d_awdens(1) ', & wdens(1), awdens(1), d_awdens(1) PRINT *, 'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ENDIF sigmaw = sigmaw + d_sigmaw wdens = wdens + d_wdens awdens = awdens + d_awdens END SUBROUTINE wake_popdyn_2 SUBROUTINE wake_popdyn_3(klon, klev, phys_sub, wk_adv, dtimesub, wgen, & wdensmin, & sigmaw, asigmaw, wdens, awdens, & !! state variables gfl, agfl, cstar, cin, wape, & rad_wk, arad_wk, irad_wk, & d_sigmaw, d_asigmaw, d_wdens, d_awdens, & !! tendencies d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, & d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd) USE lmdz_wake_ini, ONLY: wake_ini USE lmdz_wake_ini, ONLY: prt_level, RG USE lmdz_wake_ini, ONLY: stark, wdens_ref USE lmdz_wake_ini, ONLY: tau_cv, rzero, aa0 !! USE lmdz_wake_ini , ONLY: iflag_wk_pop_dyn, wdensmin USE lmdz_wake_ini, ONLY: iflag_wk_pop_dyn USE lmdz_wake_ini, ONLY: sigmad, cstart, sigmaw_max USE lmdz_wake_ini, ONLY: smallestreal IMPLICIT NONE INTEGER, INTENT(IN) :: klon, klev LOGICAL, INTENT(IN) :: phys_sub LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv REAL, INTENT(IN) :: dtimesub REAL, INTENT(IN) :: wdensmin REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes REAL, DIMENSION (klon), INTENT(INOUT) :: asigmaw !! sigma = fractional area of active wakes REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = mean wake radius REAL, DIMENSION (klon), INTENT(OUT) :: arad_wk !! r_A = wake radius of active wakes REAL, DIMENSION (klon), INTENT(OUT) :: irad_wk !! r_I = wake radius of inactive wakes REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front length per unit area REAL, DIMENSION (klon), INTENT(OUT) :: agfl !! LgA = gust front length of active wakes !! per unit area REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_asigmaw, d_wdens, d_awdens ! Some components of the tendencies of state variables REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd REAL, DIMENSION (klon), INTENT(OUT) :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_acol, d_adens_icol, d_adens_bnd !! internal variables INTEGER :: i, k REAL, DIMENSION (klon) :: iwdens, isigmaw !! inactive wake density and fractional area !! REAL, DIMENSION (klon) :: d_arad, d_irad REAL, DIMENSION (klon) :: igfl !! LgI = gust front length of inactive wakes !! per unit area REAL, DIMENSION (klon) :: s_wk !! mean area of individual wakes REAL, DIMENSION (klon) :: as_wk !! mean area of individual active wakes REAL, DIMENSION (klon) :: is_wk !! mean area of individual inactive wakes REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes REAL :: tau_wk_inv_min REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes REAL :: d_wdens_targ, d_sigmaw_targ !! Equations !! --------- !! Gust fronts: !! Lg_A = 2 pi r_A A !! Lg_I = 2 pi r_I I !! Lg = 2 pi r D !! !! Areas: !! s = pi r^2 !! s_A = pi r_A^2 !! s_I = pi r_I^2 !! !! Life expectancy: !! tau_I = 3 C* ((C*/C*t)^3/2 - 1) / r_I !! !! Time deratives: !! dD/dt = B - (D-A)/tau_I - 2 Lg C* D !! dA/dt = B - A/tau_A + 2 Lg_I C* (D-A) - 2 Lg_A C* A !! dsigma/dt = B a0 - sigma_I/tau_I + Lg C* - 2 Lg_I C* (D-A) (2 s_I - a0) !! dsigma_A/dt = B a0 - sigma_A/tau_A + Lg_A C* + (Lg_A I + Lg_I A) C* s_I + 2 Lg_I C* I a0 !! DO i = 1, klon IF (wk_adv(i)) THEN iwdens(i) = wdens(i) - awdens(i) isigmaw(i) = sigmaw(i) - asigmaw(i) arad_wk(i) = max(sqrt(asigmaw(i) / (3.14 * awdens(i))), rzero) irad_wk(i) = max(sqrt((sigmaw(i) - asigmaw(i)) / & (3.14 * max(smallestreal, (wdens(i) - awdens(i))))), rzero) rad_wk(i) = (awdens(i) * arad_wk(i) + (wdens(i) - awdens(i)) * irad_wk(i)) / wdens(i) s_wk(i) = 3.14 * rad_wk(i)**2 as_wk(i) = 3.14 * arad_wk(i)**2 is_wk(i) = 3.14 * irad_wk(i)**2 gfl(i) = 2. * sqrt(3.14 * wdens(i) * sigmaw(i)) agfl(i) = 2. * sqrt(3.14 * awdens(i) * asigmaw(i)) igfl(i) = gfl(i) - agfl(i) ENDIF ENDDO DO i = 1, klon IF (wk_adv(i)) THEN tau_wk_inv(i) = max((3. * cstar(i)) / (irad_wk(i) * ((cstar(i) / cstart)**1.5 - 1)), 0.) tau_wk_inv_min = min(tau_wk_inv(i), 1. / dtimesub) tau_prime(i) = tau_cv d_sig_gen(i) = wgen(i) * aa0 d_sig_death(i) = - isigmaw(i) * tau_wk_inv_min d_sig_col(i) = - 2. * igfl(i) * cstar(i) * iwdens(i) * (2. * is_wk(i) - aa0) d_sig_spread(i) = gfl(i) * cstar(i) d_sig_gen(i) = d_sig_gen(i) * dtimesub d_sig_death(i) = d_sig_death(i) * dtimesub d_sig_col(i) = d_sig_col(i) * dtimesub d_sig_spread(i) = d_sig_spread(i) * dtimesub d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) CALL iophys_ecrit('d_sigmaw0', 1, 'd_sigmaw0', '', d_sigmaw) END IF d_sigmaw_targ = max(d_sigmaw(i), sigmad - sigmaw(i)) !! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) !! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) d_sigmaw(i) = d_sigmaw_targ !! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) THEN CALL iophys_ecrit('tauwk_inv', 1, 'tau_wk_inv_min', '', tau_wk_inv_min) CALL iophys_ecrit('d_sigmaw', 1, 'd_sigmaw', '', d_sigmaw) CALL iophys_ecrit('d_sig_gen', 1, 'd_sig_gen', '', d_sig_gen) CALL iophys_ecrit('d_sig_death', 1, 'd_sig_death', '', d_sig_death) CALL iophys_ecrit('d_sig_col', 1, 'd_sig_col', '', d_sig_col) CALL iophys_ecrit('d_sig_spread', 1, 'd_sig_spread', '', d_sig_spread) CALL iophys_ecrit('d_sig_bnd', 1, 'd_sig_bnd', '', d_sig_bnd) ENDIF END IF d_asig_death(i) = - asigmaw(i) / tau_prime(i) d_asig_aicol(i) = (agfl(i) * iwdens(i) + igfl(i) * awdens(i)) * cstar(i) * is_wk(i) d_asig_iicol(i) = 2. * igfl(i) * cstar(i) * iwdens(i) * aa0 d_asig_spread(i) = agfl(i) * cstar(i) d_asig_death(i) = d_asig_death(i) * dtimesub d_asig_aicol(i) = d_asig_aicol(i) * dtimesub d_asig_iicol(i) = d_asig_iicol(i) * dtimesub d_asig_spread(i) = d_asig_spread(i) * dtimesub d_asigmaw(i) = d_sig_gen(i) + d_asig_death(i) + d_asig_aicol(i) + d_asig_iicol(i) + d_asig_spread(i) IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) CALL iophys_ecrit('d_asigmaw0', 1, 'd_asigmaw0', '', d_asigmaw) END IF d_sigmaw_targ = min(max(d_asigmaw(i), -asigmaw(i)), sigmaw(i) - asigmaw(i)) !! d_dens_bnd(i) = d_dens_bnd(i) + d_sigmaw_targ - d_sigmaw(i) d_asig_bnd(i) = d_sigmaw_targ - d_asigmaw(i) d_asigmaw(i) = d_sigmaw_targ IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) THEN CALL iophys_ecrit('d_asigmaw', 1, 'd_asigmaw', '', d_asigmaw) CALL iophys_ecrit('d_asig_death', 1, 'd_asig_death', '', d_asig_death) CALL iophys_ecrit('d_asig_aicol', 1, 'd_asig_aicol', '', d_asig_aicol) CALL iophys_ecrit('d_asig_iicol', 1, 'd_asig_iicol', '', d_asig_iicol) CALL iophys_ecrit('d_asig_spread', 1, 'd_asig_spread', '', d_asig_spread) CALL iophys_ecrit('d_asig_bnd', 1, 'd_asig_bnd', '', d_asig_bnd) ENDIF END IF d_dens_gen(i) = wgen(i) d_dens_death(i) = - iwdens(i) * tau_wk_inv_min d_dens_col(i) = - 2. * gfl(i) * cstar(i) * wdens(i) d_dens_gen(i) = d_dens_gen(i) * dtimesub d_dens_death(i) = d_dens_death(i) * dtimesub d_dens_col(i) = d_dens_col(i) * dtimesub d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i) !! d_wdens_targ = max(d_wdens(i), wdensmin - wdens(i)) !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) d_dens_bnd(i) = d_wdens_targ - d_wdens(i) d_wdens(i) = d_wdens_targ IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) THEN CALL iophys_ecrit('d_wdens', 1, 'd_wdens', '', d_wdens) CALL iophys_ecrit('d_dens_gen', 1, 'd_dens_gen', '', d_dens_gen) CALL iophys_ecrit('d_dens_death', 1, 'd_dens_death', '', d_dens_death) CALL iophys_ecrit('d_dens_col', 1, 'd_dens_col', '', d_dens_col) ENDIF END IF d_adens_death(i) = -awdens(i) / tau_prime(i) d_adens_icol(i) = 2. * igfl(i) * cstar(i) * iwdens(i) d_adens_acol(i) = - 2. * agfl(i) * cstar(i) * awdens(i) d_adens_death(i) = d_adens_death(i) * dtimesub d_adens_icol(i) = d_adens_icol(i) * dtimesub d_adens_acol(i) = d_adens_acol(i) * dtimesub d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i) IF (CPPKEY_IOPHYS_WK) THEN IF (phys_sub) THEN CALL iophys_ecrit('d_awdens', 1, 'd_awdens', '', d_awdens) CALL iophys_ecrit('d_adens_death', 1, 'd_adens_death', '', d_adens_death) CALL iophys_ecrit('d_adens_icol', 1, 'd_adens_icol', '', d_adens_icol) CALL iophys_ecrit('d_adens_acol', 1, 'd_adens_acol', '', d_adens_acol) ENDIF END IF d_wdens_targ = min(max(d_awdens(i), -awdens(i)), wdens(i) - awdens(i)) !! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) d_adens_bnd(i) = d_wdens_targ - d_awdens(i) d_awdens(i) = d_wdens_targ !! d_irad(i) = (d_sigmaw(i)-d_asigmaw(i)-isigmaw(i)*(d_wdens(i)-awdens(i))/iwdens(i)) / & !! max(smallestreal,(2.*3.14*iwdens(i)*irad_wk(i))) !! d_arad(i) = (d_asigmaw(i)-asigmaw(i)*d_awdens(i)/awdens(i)) / & !! max(smallestreal,(2.*3.14*awdens(i)*arad_wk(i))) !! d_irad(i) = d_irad(i)*dtimesub !! d_arad(i) = d_arad(i)*dtimesub !! CALL iophys_ecrit('d_irad',1,'d_irad','m',d_irad) !! CALL iophys_ecrit('d_airad',1,'d_arad','m',d_arad) !! ENDIF ENDDO IF (prt_level >= 10) THEN PRINT *, 'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1) ', & cstar(1), cstar(1) / cstart, rad_wk(1), tau_wk_inv(1), gfl(1) PRINT *, 'wake, wdens(1), awdens(1), d_awdens(1) ', & wdens(1), awdens(1), d_awdens(1) PRINT *, 'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ENDIF sigmaw = sigmaw + d_sigmaw asigmaw = asigmaw + d_asigmaw wdens = wdens + d_wdens awdens = awdens + d_awdens END SUBROUTINE wake_popdyn_3 END MODULE lmdz_wake