Index: LMDZ6/trunk/libf/phylmd/lmdz_wake.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_wake.f90 (revision 5760) +++ LMDZ6/trunk/libf/phylmd/lmdz_wake.f90 (revision 5761) @@ -5,5 +5,6 @@ IMPLICIT NONE; PRIVATE - LOGICAL, PARAMETER :: phys_sub=.false. +!!! LOGICAL, PARAMETER :: phys_sub=.false. + LOGICAL, PARAMETER :: phys_sub=.true. LOGICAL :: first_call=.true. !$OMP THREADPRIVATE(first_call) @@ -42,5 +43,6 @@ dth, hw, wape, fip, gfl, & dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, & - dtke, dqke, omg, dp_deltomg, wkspread, cstar, & +!!! dtke, dqke, omg, dp_deltomg, wkspread, cstar, & ! changes in notation + d_deltat_dcv, d_deltaq_dcv, omg, dp_deltomg, wkspread, cstar, & d_deltat_gw, & ! tendencies d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2) ! tendencies @@ -96,6 +98,6 @@ ! wdens : densite de poches ! omgbdth: flux of Delta_Theta transported by LS omega - ! dtKE : differential heating (wake - unpertubed) - ! dqKE : differential moistening (wake - unpertubed) + ! d_deltat_dcv : differential heating (wake - unpertubed) + ! d_deltat_dcv : differential moistening (wake - unpertubed) ! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) ! dp_deltomg : vertical gradient of omg (s-1) @@ -168,5 +170,5 @@ ! -------------------- - INTEGER, INTENT(IN) :: klon,klev + INTEGER, INTENT(IN) :: klon,klev INTEGER, DIMENSION (klon), INTENT(IN) :: znatsurf REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi @@ -197,5 +199,6 @@ 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) :: dtke, dqke + REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_dcv, d_deltaq_dcv REAL, DIMENSION (klon, klev), INTENT(OUT) :: wkspread ! unused (jyg) REAL, DIMENSION (klon, klev), INTENT(OUT) :: omgbdth, omg @@ -267,5 +270,13 @@ ! Sub-timestep tendencies and related variables REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw + REAL, DIMENSION (klon, klev) :: d_deltat_dadv, d_deltaq_dadv + REAL, DIMENSION (klon, klev) :: d_deltat_lsadv, d_deltaq_lsadv + REAL, DIMENSION (klon, klev) :: d_deltat_entrp + REAL, DIMENSION (klon, klev) :: d_deltat_spread, d_deltaq_spread + REAL, DIMENSION (klon, klev) :: d_tb, d_qb + REAL, DIMENSION (klon, klev) :: d_tb_dadv, d_qb_dadv + REAL, DIMENSION (klon, klev) :: d_tb_spread, d_qb_spread + 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 @@ -279,4 +290,7 @@ LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw + ! Tests + REAL, DIMENSION (klon, klev) :: c5p + ! Autres variables internes INTEGER ::isubstep, k, i, igout @@ -288,4 +302,6 @@ REAL :: d_sigmaw_targ REAL :: d_wdens_targ + + REAL, DIMENSION (klon) :: dsigspread !rate of change of sigmaw due to spreading REAL, DIMENSION (klon) :: sum_thx, sum_tx, sum_qx, sum_thvx @@ -317,5 +333,4 @@ REAL, DIMENSION (klon, klev) :: omgbdq - REAL, DIMENSION (klon) :: ff, gg REAL, DIMENSION (klon) :: wape2, cstar2, heff @@ -327,6 +342,7 @@ REAL, DIMENSION (klon) :: death_rate !! REAL, DIMENSION (klon) :: nat_rate - REAL, DIMENSION (klon, klev) :: entr - REAL, DIMENSION (klon, klev) :: detr + REAL, DIMENSION (klon, klev) :: entr ! total entrainment into wakes (spread + birth) + REAL, DIMENSION (klon, klev) :: entr_p ! entrainment into wakes (due to births) + REAL, DIMENSION (klon, klev) :: detr ! detrainment from wakes (due to deaths) REAL, DIMENSION(klon) :: sigmaw_in, asigmaw_in ! pour les prints @@ -341,4 +357,8 @@ REAL, DIMENSION (klon) :: omega REAL, DIMENSION (klon) :: h_zzz + + !! Bidouilles + REAL :: iwkadv + REAL :: iokqxqw !print*,'WAKE LJYFz' @@ -438,8 +458,35 @@ dqls(:,:) = 0. d_deltat_gw(:,:) = 0. + + detr(:,:) = 0. + entr(:,:) = 0. + entr_p(:,:) = 0. + c5p(:,:) = 0. + + th1(:,:) = 0. + th2(:,:) = 0. + q1(:,:) = 0. + q2(:,:) = 0. + d_tb(:,:) = 0. + d_tb_dadv(:,:) = 0. + d_tb_spread(:,:) = 0. + d_qb(:,:) = 0. + d_qb_dadv(:,:) = 0. + d_qb_spread(:,:) = 0. + d_deltatw(:,:) = 0. + d_deltat_dadv (:,:) = 0. + d_deltat_lsadv (:,:) = 0. + d_deltat_dcv (:,:) = 0. + d_deltat_spread(:,:) = 0. + d_deltaqw(:,:) = 0. + d_deltaq_dadv(:,:) = 0. + d_deltaq_lsadv(:,:) = 0. + d_deltaq_dcv(:,:) = 0. + d_deltaq_spread(:,:) = 0. + d_deltatw2(:,:) = 0. d_deltaqw2(:,:) = 0. @@ -563,6 +610,6 @@ tx(:,:) = 0. qx(:,:) = 0. -dtke(:,:) = 0. -dqke(:,:) = 0. +d_deltat_dcv(:,:) = 0. +d_deltaq_dcv(:,:) = 0. wkspread(:,:) = 0. omgbdth(:,:) = 0. @@ -883,4 +930,5 @@ ! ! ------------------------------------------------------------------------ + ! ! wk_adv is the logical flag enabling wake evolution in the time advance ! loop @@ -888,4 +936,15 @@ wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1. END DO +IF (CPPKEY_IOPHYS_WK) THEN + IF (phys_sub) THEN + iwkadv=0. + IF (wk_adv(1)) iwkadv=1. + iokqxqw=0. + IF (ok_qx_qw(1)) iokqxqw=1. + CALL iophys_ecrit('iwkadv',1,'iwkadv','',iwkadv) + CALL iophys_ecrit('iokqxqw',1,'iokqxqw','',iokqxqw) + CALL iophys_ecrit('alpha',1,'alpha','',alpha(1)) + ENDIF +END IF IF (prt_level>=10) THEN PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', & @@ -1021,11 +1080,10 @@ IF (phys_sub) THEN CALL iophys_ecrit('ptop',1,'ptop','Pa',ptop) + CALL iophys_ecrit('wape',1,'wape','J/kg',wape) + CALL iophys_ecrit('wgen',1,'wgen','1/(s.m2)',wgen) 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 @@ -1041,4 +1099,11 @@ d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd ) +IF (CPPKEY_IOPHYS_WK) THEN + IF (phys_sub) THEN + 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 sigmaw=sigmaw-d_sigmaw asigmaw=asigmaw-d_asigmaw @@ -1313,4 +1378,26 @@ ENDIF + + ! ----------------------------------------------------------------- + ! Compute redistribution (advective) term + +! rate of change of sigmaw due to spreading + dsigspread(:) = gfl(:)*cstar(:) + + CALL wake_dadv(klon, klev, dtimesub, ph, ppi, wk_adv, kupper, & + omg, dp_deltomg, sigmaw, dsigspread, & + th2, th1, q2, q1, & + d_deltat_dadv, d_deltaq_dadv, d_tb_dadv, d_qb_dadv) + + ! For the difference fields: convert to change per second in order to combine with the + ! other terms (d_deltat_ls, d_deltat_cv, d_deltat_gw) + d_deltat_dadv(:,:) = d_deltat_dadv(:,:)/dtimesub + d_deltaq_dadv(:,:) = d_deltaq_dadv(:,:)/dtimesub +! + ! For the mean fields: tb and qb the computation of the tendencies due to wakes is + ! already complete. + d_tb(:,:) = d_tb_dadv(:,:) + d_qb(:,:) = d_qb_dadv(:,:) + ! ----------------------------------------------------------------- DO k = 1, klev-1 @@ -1318,43 +1405,11 @@ 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)- & + d_deltat_lsadv(i, k) = 1./(ph(i,k)-ph(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)- & + + d_deltaq_lsadv(i, k) = 1./(ph(i,k)-ph(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 @@ -1364,7 +1419,10 @@ ! ------------------------------------------------------------------ - 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) +!! IF (prt_level>=10) THEN + IF (prt_level>=10 .and. wk_adv(igout)) THEN + PRINT *, 'wake-4.3, d_deltat_dadv(igout,k) ', (k,d_deltat_dadv(igout,k), k=1,klev) + PRINT *, 'wake-4.3, d_deltat_lsadv(igout,k) ', (k,d_deltat_lsadv(igout,k), k=1,klev) + PRINT *, 'wake-4.3, d_deltaq_dadv(igout,k) ', (k,d_deltaq_dadv(igout,k), k=1,klev) + PRINT *, 'wake-4.3, d_deltaq_lsadv(igout,k) ', (k,d_deltaq_lsadv(igout,k), k=1,klev) ENDIF @@ -1376,5 +1434,5 @@ 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) + entr_p(i,k) = d_sig_gen(i) ENDIF ENDDO @@ -1386,5 +1444,5 @@ detr(i, k) = 0. - entr(i, k) = 0. + entr_p(i, k) = 0. ENDIF ENDDO @@ -1426,6 +1484,9 @@ ! . /(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))) + !! Differential heating (d_deltat_dcv) and moistening (d_deltaq_dcv) by deep convection + d_deltat_dcv(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i))) +!! dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i))) ! supprime + d_deltaq_dcv(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i))) +!! dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i))) ! supprime ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k) @@ -1440,5 +1501,5 @@ !! sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) !! - entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + & + entr(i, k) = entr_p(i,k) + gfl(i)*cstar(i) + & sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) !>jyg @@ -1458,5 +1519,4 @@ dtimesub ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k) - ff(i) = d_deltatw(i, k)/dtimesub ! Sans GW @@ -1471,13 +1531,18 @@ ! GW formule 2 + !! Entrainment due to spread is supposed to be included in the differential advection + !! term (d_deltat_dadv); hence only the entrainment due to population dynamics (entr_p) + !! appears in the expression of d_deltatw. 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) - & +!!! d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - & ! nouvelle notation + d_deltatw(i, k) = dtimesub*(d_deltat_dadv(i,k)+d_deltat_lsadv(i,k)+d_deltat_dcv(i,k) - & + entr_p(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) - & +!!! (ff(i)+dtke(i,k) - & ! nouvelle notation + (d_deltat_dadv(i,k)+d_deltat_lsadv(i,k)+d_deltat_dcv(i,k) - & + entr_p(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) ) @@ -1485,9 +1550,11 @@ 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) - & + c5p(i,k) = - entr_p(i,k)*deltatw(i,k)/sigmaw(i) + + !! Entrainment due to spread is supposed to be included in the differential advection + !! term (d_deltaq_dadv); hence only the entrainment due to population dynamics (entr_p) + !! appears in the expression of d_deltaqw. + d_deltaqw(i, k) = dtimesub*(d_deltaq_dadv(i,k)+d_deltaq_lsadv(i,k)+d_deltaq_dcv(i,k) - & + entr_p(i,k)*deltaqw(i,k)/sigmaw(i) - & (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i))) ! cc @@ -1501,4 +1568,32 @@ END DO +!! IF (prt_level>=10) THEN + IF (prt_level>=10 .and. wk_adv(igout)) THEN + PRINT *, 'wake-4.4, isubstep= ', isubstep,' c5p(igout,k) ', (k,c5p(igout,k), k=1,klev) + PRINT *, 'wake-4.4, isubstep= ', isubstep,' deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev) + PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltat_dcv(igout,k) ', (k,d_deltat_dcv(igout,k), k=1,klev) + PRINT *, 'wake-4.4, isubstep= ', isubstep,' tgw(igout,k)*deltatw(igout,k) ', (k,tgw(igout,k)*deltatw(igout,k), k=1,klev) + PRINT *, 'wake-4.4, isubstep= ', isubstep,' death_rate(igout) ', death_rate(igout) + PRINT *, 'wake-4.4, isubstep= ', isubstep,' detr(igout,k) ', (k,detr(igout,k), k=1,klev) + PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev) + PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev) + PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_tb(igout,k) ', (k,d_tb(igout,k), k=1,klev) + PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_qb(igout,k) ', (k,d_qb(igout,k), k=1,klev) + ENDIF + +! +IF (CPPKEY_IOPHYS_WK) THEN + IF (phys_sub) THEN + DO k = 1,klev + d_deltat_entrp(:,k) = - entr_p(:,k)*deltatw(:,k)/sigmaw(:) + ENDDO + CALL iophys_ecrit('d_deltatw',klev,'d_deltatw','K/s',d_deltatw(:,1:klev)) + CALL iophys_ecrit('d_deltat_dadv',klev,'d_deltat_dadv','K/s',d_deltat_dadv(:,1:klev)) + CALL iophys_ecrit('d_deltat_lsadv',klev,'d_deltat_lsadv','K/s',d_deltat_lsadv(:,1:klev)) + CALL iophys_ecrit('d_deltat_dcv',klev,'d_deltat_dcv','K/s',d_deltat_dcv(:,1:klev)) + CALL iophys_ecrit('d_deltat_entrp',klev,'d_deltat_entrp','K/s',d_deltat_entrp(:,1:klev)) + + ENDIF +END IF ! Scale tendencies so that water vapour remains positive in w and x. @@ -1822,4 +1917,6 @@ 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) +!!print *,'XXXXwake wape(i), hw0(i), av_dth(i), av_thx(i), av_dq(i), av_qx(i), av_thvx(i) ', & +!! wape(i), hw0(i), av_dth(i), av_thx(i), av_dq(i), av_qx(i), av_thvx(i) END IF END DO @@ -1879,4 +1976,5 @@ IF (CPPKEY_IOPHYS_WK) THEN IF (.not.phys_sub) CALL iophys_ecrit('wape_b',1,'wape_b','J/kg',wape) + IF (.not.phys_sub) CALL iophys_ecrit('alpha_mod',1,'alpha_modulation','-',alpha_tot) END IF IF (prt_level>=10) THEN @@ -2136,16 +2234,4 @@ 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) @@ -2281,5 +2367,20 @@ IF (CPPKEY_IOPHYS_WK) THEN IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape) -END IF + IF (iflag_wk_pop_dyn >= 3) 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) + ENDIF ! (.not.phys_sub) + ENDIF ! (iflag_wk_pop_dyn >= 3) +END IF !(CPPKEY_IOPHYS_WK) @@ -2363,11 +2464,11 @@ ! 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 + INTEGER, INTENT(IN) :: nl, nlon + REAL, DIMENSION(nlon, nl), INTENT(IN) :: qb, d_qb + REAL, DIMENSION(nlon, nl), INTENT(IN) :: deltaqw, d_deltaqw + REAL, DIMENSION(nlon), INTENT(IN) :: sigmaw, d_sigmaw + LOGICAL, DIMENSION(nlon), INTENT(IN) :: wk_adv ! Output - REAL alpha(nlon) + REAL, DIMENSION(nlon), INTENT(INOUT) :: alpha ! Internal variables REAL zeta(nlon, nl) @@ -3202,5 +3303,4 @@ 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 @@ -3229,4 +3329,38 @@ !! +! Initialization + tau_wk_inv(:) = 0. +! Initialization of output variables + rad_wk(:) = 0. + arad_wk(:) = 0. + irad_wk(:) = 0. + gfl(:) = 0. + agfl(:) = 0. +! + d_wdens(:) = 0. + d_awdens(:) = 0. + d_sigmaw(:) = 0. + d_asigmaw (:) = 0. +! + d_sig_gen(:) = 0. + d_sig_death(:) = 0. + d_sig_col(:) = 0. + d_sig_spread(:) = 0. + d_sig_bnd(:) = 0. + d_asig_death(:) = 0. + d_asig_aicol(:) = 0. + d_asig_iicol(:) = 0. + d_asig_spread(:) = 0. + d_asig_bnd(:) = 0. + d_dens_gen(:) = 0. + d_dens_death(:) = 0. + d_dens_col(:) = 0. + d_dens_bnd(:) = 0. + d_adens_death(:) = 0. + d_adens_icol(:) = 0. + d_adens_acol(:) = 0. + d_adens_bnd(:) = 0. + + DO i = 1, klon IF (wk_adv(i)) THEN @@ -3252,10 +3386,14 @@ 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) +! print *,'ZZZZpopdyn3 wgen(1) ',wgen(1) +! print *,'ZZZZpopdyn3 cstar(1) ',cstar(1) +! print *,'ZZZZpopdyn3 isigmaw(1) ',isigmaw(1) +! print *,'ZZZZpopdyn3 gfl(1) ',gfl(1) +!! tau_wk_inv(i) = max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.) + tau_wk_inv(i) = min(max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.), 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_death(i) = - isigmaw(i)*tau_wk_inv(i) d_sig_col(i) = - 2.*igfl(i)*cstar(i)*iwdens(i)*(2.*is_wk(i)-aa0) d_sig_spread(i) = gfl(i)*cstar(i) @@ -3266,8 +3404,4 @@ 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)) @@ -3277,15 +3411,4 @@ 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) @@ -3298,24 +3421,11 @@ 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_death(i) = - iwdens(i)*tau_wk_inv(i) d_dens_col(i) = - 2.*gfl(i)*cstar(i)*wdens(i) ! @@ -3329,13 +3439,4 @@ 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) @@ -3346,12 +3447,4 @@ 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) @@ -3370,4 +3463,39 @@ ENDIF ENDDO +IF (CPPKEY_IOPHYS_WK) THEN + IF (phys_sub) THEN + call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw) +! + call iophys_ecrit('cstar',1,'cstar','',cstar) + call iophys_ecrit('wgen_pd3',1,'wgen_popdyn3','',wgen) + call iophys_ecrit('tauwk_inv',1,'tau_wk_inv','',tau_wk_inv) + 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) +! + call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw) +! + 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) +! + 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) +! + 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 + IF (prt_level >= 10) THEN @@ -3386,4 +3514,458 @@ RETURN END SUBROUTINE wake_popdyn_3 + + SUBROUTINE wake_dadv(klon, klev, dtime, ph, ppi, wk_adv, kupper, & + deltomg, dp_deltomg, sigmaw, dsigspread, & + thw, thx, qw, qx, & + d_deltat_dadv, d_deltaq_dadv, d_tb_dadv, d_qb_dadv) + + USE lmdz_wake_ini , ONLY : flag_dadv_implicit + +IMPLICIT NONE + + INTEGER, INTENT(IN) :: klon, klev + REAL, INTENT(IN) :: dtime + REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph + REAL, DIMENSION (klon, klev), INTENT(IN) :: ppi + LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv + INTEGER, DIMENSION (klon), INTENT(IN) :: kupper + REAL, DIMENSION (klon, klev), INTENT(IN) :: deltomg + REAL, DIMENSION (klon, klev), INTENT(IN) :: dp_deltomg + REAL, DIMENSION (klon), INTENT(IN) :: sigmaw + REAL, DIMENSION (klon), INTENT(IN) :: dsigspread + REAL, DIMENSION (klon, klev), INTENT(IN) :: thw ! component # 1 + REAL, DIMENSION (klon, klev), INTENT(IN) :: thx ! component # 2 + REAL, DIMENSION (klon, klev), INTENT(IN) :: qw ! component # 1 + REAL, DIMENSION (klon, klev), INTENT(IN) :: qx ! component # 2 + + REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_dadv + REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltaq_dadv + REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_tb_dadv + REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_qb_dadv + +! Internal variables + INTEGER :: i, k + REAL, DIMENSION (klon, klev) :: entr_s ! entrainment into wakes due to spread + REAL, DIMENSION (klon, klev) :: thb, qb + REAL, DIMENSION (klon, klev) :: delta_th, delta_q + +! Tests + +! Arrays used in the implicit scheme + REAL, DIMENSION (klon) :: rr11, rr12, rr21, rr22 + + REAL, DIMENSION (klon, klev) :: aa11, aa12, aa21, aa22 + REAL, DIMENSION (klon, klev) :: bb11, bb12, bb21, bb22 + REAL, DIMENSION (klon, klev) :: cc11, cc12, cc21, cc22 + + REAL, DIMENSION (klon, klev) :: alpha11, alpha12, alpha21, alpha22 + REAL, DIMENSION (klon, klev) :: beta11, beta12, beta21, beta22 + REAL, DIMENSION (klon, klev) :: gamma11, gamma12, gamma21, gamma22 + REAL, DIMENSION (klon, klev) :: ai11, ai12, ai21, ai22 ! inverse of alpha + + REAL, DIMENSION (klon, klev) :: xt1, xt2, xq1, xq2 + REAL, DIMENSION (klon, klev) :: yt1, yt2, yq1, yq2 + REAL, DIMENSION (klon, klev) :: zt1, zt2, zq1, zq2 + REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2 + + REAL :: coef, det + + REAL, DIMENSION (klon,klev) :: xt1inv, xt2inv, xq1inv, xq2inv + +! Arrays used in the explicit scheme (vertical gradients) + REAL, DIMENSION (klon, klev) :: d_thx, d_qx + REAL, DIMENSION (klon, klev) :: d_thw, d_qw + REAL, DIMENSION (klon, klev) :: d_dth, d_dq + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + entr_s(:,:) = 0. + delta_th(:,:) = 0. + + d_deltat_dadv(:,:) = 0. + d_deltaq_dadv(:,:) = 0. + d_tb_dadv(:,:) = 0. + d_qb_dadv(:,:) = 0. + + + rr11(:) = sigmaw(:) + rr12(:) = 1.-sigmaw(:) + rr21(:) = 1. + rr22(:) = -1. + + DO k = 1, klev + DO i = 1,klon + IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN + thb(i,k) = rr11(i)*thw(i,k)+rr12(i)*thx(i,k) + delta_th(i,k) = rr21(i)*thw(i,k)+rr22(i)*thx(i,k) + + qb(i,k) = rr11(i)*qw(i,k)+rr12(i)*qx(i,k) + delta_q(i,k) = rr21(i)*qw(i,k)+rr22(i)*qx(i,k) + ENDIF + ENDDO + ENDDO + + DO i = 1, klon + entr_s(i,klev) = 0. + ENDDO + + DO k = 1, klev-1 + DO i = 1,klon + IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN +!! entr_s(i,k) = dsigspread(i) - sigmaw(i)*(1.-sigmaw(i))*(deltomg(i,k+1)-deltomg(i,k)) / & +!! (ph(i,k)-ph(i,k+1)) + + entr_s(i,k) = dsigspread(i) + sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i,k) +!! print *,'dadv, k, dp_deltomg(i,k), (deltomg(i,k)-deltomg(i,k+1))/(ph(i,k)-ph(i,k+1)) ', & +!! k, dp_deltomg(i,k), (deltomg(i,k)-deltomg(i,k+1))/(ph(i,k)-ph(i,k+1)) + + ENDIF + ENDDO + ENDDO + +! ------------------------------------------------------------------------------------------- +! Depending on flag_dadv_implicit, use implicit upstream scheme or explicit upstream scheme +! ------------------------------------------------------------------------------------------- + + IF (flag_dadv_implicit) THEN + +! Implicit scheme : solve for d_deltat_dadv and d_tb_dadv +! (and similarly for d_deltaq_dadv and d_qb_dadv). +! The system to be inverted is block-tridiagonal with 2x2 blocks. +! ----------------------------------------------------------------------------------------- + +! Matrix indexing: Theta_w Theta_x +! +! / +! Theta_b | A11 A12 | +! | | +! delta_theta | A21 A22 | +! / +! Tridiagonal matrix +! / +! | aa(1) bb(1) 0 | +! | cc(2) aa(2) bb(2) 0 | +! | 0 cc(3) aa(3) bb(3) | +! | | +! . . . . +! +! . . . . +! | | +! | cc(n-2) aa(n-2) bb(n-2) 0 | +! | 0 cc(n-1) aa(n-1) bb(n-1) | +! 0 cc(n) aa(n) / +! ----------------------------------------------------------------------------------------- + +!! Building the tridiagonal matrix + DO i = 1,klon + IF (wk_adv(i)) THEN + k = kupper(i) + coef = dtime/(ph(i,k)-ph(i,k+1)) + aa11(i,k) = rr11(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k) + aa12(i,k) = rr12(i) + aa21(i,k) = rr21(i)+coef*( dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) + & + (1.-sigmaw(i))*deltomg(i,k) ) + aa22(i,k) = rr22(i)+coef*(-dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) - & + deltomg(i,k) ) + + cc11(i,k) = 0. + cc12(i,k) = -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k) + cc21(i,k) = 0. + cc22(i,k) = coef*sigmaw(i)*deltomg(i,k) + ENDIF ! (wk_adv(i)) + ENDDO + DO k = 2, klev-1 + DO i = 1,klon + IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN + coef = dtime/(ph(i,k)-ph(i,k+1)) + aa11(i,k) = rr11(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k) + aa12(i,k) = rr12(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k+1) + aa21(i,k) = rr21(i)+coef*( dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) + & + (1.-sigmaw(i))*deltomg(i,k) ) + aa22(i,k) = rr22(i)+coef*(-dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) + & + (1.-sigmaw(i))*(deltomg(i,k+1)-deltomg(i,k)) - & + sigmaw(i)*deltomg(i,k) ) + + bb11(i,k) = -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k+1) + bb12(i,k) = 0. + bb21(i,k) = -coef*(1.-sigmaw(i))*deltomg(i,k+1) + bb22(i,k) = 0. + + cc11(i,k) = 0. + cc12(i,k) = -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k) + cc21(i,k) = 0. + cc22(i,k) = coef*sigmaw(i)*deltomg(i,k) + ENDIF ! (wk_adv(i) .AND. k<=kupper(i)) + ENDDO + ENDDO + DO i = 1,klon + IF (wk_adv(i)) THEN + coef = dtime/(ph(i,1)-ph(i,2)) + aa11(i,1) = rr11(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,1) + aa12(i,1) = rr12(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,2) + aa21(i,1) = rr21(i)+coef*( dsigspread(i)/sigmaw(i)*(ph(i,1)-ph(i,2)) + & + (1.-sigmaw(i))*deltomg(i,1) ) + aa22(i,1) = rr22(i)+coef*(-dsigspread(i)/sigmaw(i)*(ph(i,1)-ph(i,2)) + & + (1.-sigmaw(i))*(deltomg(i,2)-deltomg(i,1)) - & + sigmaw(i)*deltomg(i,1) ) + + bb11(i,1) = -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,2) + bb12(i,1) = 0. + bb21(i,1) = -coef*(1.-sigmaw(i))*deltomg(i,2) + bb22(i,1) = 0. + ENDIF ! (wk_adv(i)) + ENDDO + +!! printing the tridiagonal matrix +!!! First row +!! k = 1 +!! print 1789, k, aa11(1,1), aa12(1,1), bb11(1,1), bb12(1,1) +!! print 1789, k, aa21(1,1), aa22(1,1), bb21(1,1), bb22(1,1) +!!1789 FORMAT(1X, I3, 3(4X, 2E13.5)) +!! coef = dtime/(ph(1,k)-ph(1,k+1)) +!! print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,1), deltomg(1,2) ', & +!! rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,1), deltomg(1,2) +!! +!!! Rows 2 to klev-1 +!! DO k = 2, klev-1 +!! print 1789, k, cc11(1,k), cc12(1,k), aa11(1,k), aa12(1,k), bb11(1,k), bb12(1,k) +!! print 1789, k, cc21(1,k), cc22(1,k), aa21(1,k), aa22(1,k), bb21(1,k), bb22(1,k) +!! coef = dtime/(ph(1,k)-ph(1,k+1)) +!! print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,k), deltomg(1,k+1) ', & +!! rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,k), deltomg(1,k+1) +!! ENDDO +!! +!!! Row klev +!! print 1789, klev, cc11(1,klev), cc12(1,klev), aa11(1,klev), aa12(1,klev) +!! print 1789, klev, cc21(1,klev), cc22(1,klev), aa21(1,klev), aa22(1,klev) +!! coef = dtime/(ph(1,klev)-ph(1,klev+1)) +!! print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,klev) ', & +!! rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,klev) + + +!! Downward loop + + xt1(:,:) = thb(:,:) + xt2(:,:) = delta_th(:,:) + xq1(:,:) = qb(:,:) + xq2(:,:) = delta_q(:,:) + + DO i = 1,klon + IF (wk_adv(i)) THEN + k = kupper(i) + alpha11(:,k)=aa11(:,k) + alpha12(:,k)=aa12(:,k) + alpha21(:,k)=aa21(:,k) + alpha22(:,k)=aa22(:,k) + beta11(:,k)=0. + beta12(:,k)=0. + beta21(:,k)=0. + beta22(:,k)=0. + yt1(i,k) = xt1(i,k) + yt2(i,k) = xt2(i,k) + yq1(i,k) = xq1(i,k) + yq2(i,k) = xq2(i,k) + ENDIF ! (wk_adv(i)) + ENDDO + DO i = 1,klon + IF (wk_adv(i)) THEN + k = kupper(i) + det=alpha11(i,k)*alpha22(i,k) - alpha12(i,k)*alpha21(i,k) + ai11(i,k)= alpha22(i,k)/det + ai12(i,k)=-alpha12(i,k)/det + ai21(i,k)=-alpha21(i,k)/det + ai22(i,k)= alpha11(i,k)/det + zt1(i,k) = ai11(i,k)*yt1(i,k) + ai12(i,k)*yt2(i,k) + zt2(i,k) = ai21(i,k)*yt1(i,k) + ai22(i,k)*yt2(i,k) + zq1(i,k) = ai11(i,k)*yq1(i,k) + ai12(i,k)*yq2(i,k) + zq2(i,k) = ai21(i,k)*yq1(i,k) + ai22(i,k)*yq2(i,k) + ENDIF ! (wk_adv(i)) + ENDDO + + DO k = klev, 2, -1 + DO i = 1,klon + IF (wk_adv(i) .AND. k<=kupper(i)) THEN + gamma11(i,k) = ai11(i,k)*cc11(i,k) + ai12(i,k)*cc21(i,k) + gamma12(i,k) = ai11(i,k)*cc12(i,k) + ai12(i,k)*cc22(i,k) + gamma21(i,k) = ai21(i,k)*cc11(i,k) + ai22(i,k)*cc21(i,k) + gamma22(i,k) = ai21(i,k)*cc12(i,k) + ai22(i,k)*cc22(i,k) + + alpha11(i,k-1) = aa11(i,k-1) - ( bb11(i,k-1)*gamma11(i,k)+bb12(i,k-1)*gamma21(i,k) ) + alpha12(i,k-1) = aa12(i,k-1) - ( bb11(i,k-1)*gamma12(i,k)+bb12(i,k-1)*gamma22(i,k) ) + alpha21(i,k-1) = aa21(i,k-1) - ( bb21(i,k-1)*gamma11(i,k)+bb22(i,k-1)*gamma21(i,k) ) + alpha22(i,k-1) = aa22(i,k-1) - ( bb21(i,k-1)*gamma12(i,k)+bb22(i,k-1)*gamma22(i,k) ) + + beta11(i,k-1) = bb11(i,k-1)*ai11(i,k)+bb12(i,k-1)*ai21(i,k) + beta12(i,k-1) = bb11(i,k-1)*ai12(i,k)+bb12(i,k-1)*ai22(i,k) + beta21(i,k-1) = bb21(i,k-1)*ai11(i,k)+bb22(i,k-1)*ai21(i,k) + beta22(i,k-1) = bb21(i,k-1)*ai12(i,k)+bb22(i,k-1)*ai22(i,k) + + yt1(i,k-1) = xt1(i,k-1) - ( beta11(i,k-1)*yt1(i,k) +beta12(i,k-1)*yt2(i,k) ) + yt2(i,k-1) = xt2(i,k-1) - ( beta21(i,k-1)*yt1(i,k) +beta22(i,k-1)*yt2(i,k) ) + yq1(i,k-1) = xq1(i,k-1) - ( beta11(i,k-1)*yq1(i,k) +beta12(i,k-1)*yq2(i,k) ) + yq2(i,k-1) = xq2(i,k-1) - ( beta21(i,k-1)*yq1(i,k) +beta22(i,k-1)*yq2(i,k) ) + + det=alpha11(i,k-1)*alpha22(i,k-1) - alpha12(i,k-1)*alpha21(i,k-1) + ai11(i,k-1)= alpha22(i,k-1)/det + ai12(i,k-1)=-alpha12(i,k-1)/det + ai21(i,k-1)=-alpha21(i,k-1)/det + ai22(i,k-1)= alpha11(i,k-1)/det + + zt1(i,k-1) = ai11(i,k-1)*yt1(i,k-1)+ai12(i,k-1)*yt2(i,k-1) + zt2(i,k-1) = ai21(i,k-1)*yt1(i,k-1)+ai22(i,k-1)*yt2(i,k-1) + zq1(i,k-1) = ai11(i,k-1)*yq1(i,k-1)+ai12(i,k-1)*yq2(i,k-1) + zq2(i,k-1) = ai21(i,k-1)*yq1(i,k-1)+ai22(i,k-1)*yq2(i,k-1) + ENDIF ! (wk_adv(i) .AND. k<=kupper(i)) + ENDDO + ENDDO + +!! Upward loop + + DO i = 1,klon + IF (wk_adv(i)) THEN + th1(i,1) = zt1(i,1) + th2(i,1) = zt2(i,1) + q1(i,1) = zq1(i,1) + q2(i,1) = zq2(i,1) + + d_tb_dadv(i,1) = ( rr11(i)*(th1(i,1)-thw(i,1))+rr12(i)*(th2(i,1)-thx(i,1)) )*ppi(i,1) + d_deltat_dadv(i,1) = ( rr21(i)*(th1(i,1)-thw(i,1))+rr22(i)*(th2(i,1)-thx(i,1)) )*ppi(i,1) + d_qb_dadv(i,1) = rr11(i)*(q1(i,1) -qw(i,1)) +rr12(i)*(q2(i,1)-qx(i,1)) + d_deltaq_dadv(i,1) = rr21(i)*(q1(i,1) -qw(i,1)) +rr22(i)*(q2(i,1)-qx(i,1)) + ENDIF ! (wk_adv(i)) + ENDDO + + DO k = 2, klev + DO i = 1,klon + IF (wk_adv(i) .AND. k<=kupper(i)) THEN + th1(i,k) = zt1(i,k) - ( gamma11(i,k)*th1(i,k-1)+gamma12(i,k)*th2(i,k-1) ) + th2(i,k) = zt2(i,k) - ( gamma21(i,k)*th1(i,k-1)+gamma22(i,k)*th2(i,k-1) ) + q1(i,k) = zq1(i,k) - ( gamma11(i,k)*q1(i,k-1) +gamma12(i,k)*q2(i,k-1) ) + q2(i,k) = zq2(i,k) - ( gamma21(i,k)*q1(i,k-1) +gamma22(i,k)*q2(i,k-1) ) + + d_tb_dadv(i,k) = ( rr11(i)*(th1(i,k)-thw(i,k))+rr12(i)*(th2(i,k)-thx(i,k)) )*ppi(i,k) + d_deltat_dadv(i,k) = ( rr21(i)*(th1(i,k)-thw(i,k))+rr22(i)*(th2(i,k)-thx(i,k)) )*ppi(i,k) + d_qb_dadv(i,k) = rr11(i)*(q1(i,k)-qw(i,k)) +rr12(i)*(q2(i,k)-qx(i,k)) + d_deltaq_dadv(i,k) = rr21(i)*(q1(i,k)-qw(i,k)) +rr22(i)*(q2(i,k)-qx(i,k)) + ENDIF ! (wk_adv(i) .AND. k<=kupper(i)) + ENDDO + ENDDO + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!! Verification de l'inversion !!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! DO i = 1,klon +!! xt1inv(i,1) = aa11(i,1)*th1(i,1) + aa12(i,1)*th2(i,1) + bb11(i,1)*th1(i,2) + bb12(i,1)*th2(i,2) +!! xt2inv(i,1) = aa21(i,1)*th1(i,1) + aa22(i,1)*th2(i,1) + bb21(i,1)*th1(i,2) + bb22(i,1)*th2(i,2) +!! xq1inv(i,1) = aa11(i,1)*q1(i,1) + aa12(i,1)*q2(i,1) + bb11(i,1)*q1(i,2) + bb12(i,1)*q2(i,2) +!! xq2inv(i,1) = aa21(i,1)*q1(i,1) + aa22(i,1)*q2(i,1) + bb21(i,1)*q1(i,2) + bb22(i,1)*q2(i,2) +!! ENDDO +!! +!! DO k = 2, klev-1 +!! DO i = 1,klon +!! xt1inv(i,k) = aa11(i,k)*th1(i,k) + aa12(i,k)*th2(i,k) + bb11(i,k)*th1(i,k+1) + bb12(i,k)*th2(i,k+1) & +!! + cc11(i,k)*th1(i,k-1) + cc12(i,k)*th2(i,k-1) +!! xt2inv(i,k) = aa21(i,k)*th1(i,k) + aa22(i,k)*th2(i,k) + bb21(i,k)*th1(i,k+1) + bb22(i,k)*th2(i,k+1) & +!! + cc21(i,k)*th1(i,k-1) + cc22(i,k)*th2(i,k-1) +!! xq1inv(i,k) = aa11(i,k)*q1(i,k) + aa12(i,k)*q2(i,k) + bb11(i,k)*q1(i,k+1) + bb12(i,k)*q2(i,k+1) & +!! + cc11(i,k)*q1(i,k-1) + cc12(i,k)*q2(i,k-1) +!! xq2inv(i,k) = aa21(i,k)*q1(i,k) + aa22(i,k)*q2(i,k) + bb21(i,k)*q1(i,k+1) + bb22(i,k)*q2(i,k+1) & +!! + cc21(i,k)*q1(i,k-1) + cc22(i,k)*q2(i,k-1) +!! ENDDO +!! ENDDO +!! +!! DO i = 1,klon +!! xt1inv(i,klev) = aa11(i,klev)*th1(i,klev) + aa12(i,klev)*th2(i,klev) + cc11(i,klev)*th1(i,klev-1) + cc12(i,klev)*th2(i,klev-1) +!! xt2inv(i,klev) = aa21(i,klev)*th1(i,klev) + aa22(i,klev)*th2(i,klev) + cc21(i,klev)*th1(i,klev-1) + cc22(i,klev)*th2(i,klev-1) +!! xq1inv(i,klev) = aa11(i,klev)*q1(i,klev) + aa12(i,klev)*q2(i,klev) + cc11(i,klev)*q1(i,klev-1) + cc12(i,klev)*q2(i,klev-1) +!! xq2inv(i,klev) = aa21(i,klev)*q1(i,klev) + aa22(i,klev)*q2(i,klev) + cc21(i,klev)*q1(i,klev-1) + cc22(i,klev)*q2(i,klev-1) +!! ENDDO +!! +!! DO k = 1, 20 +!! IF (abs(xt1inv(1,k)-xt1(1,k)) .GT. 1.e-15*xt1(1,k) ) THEN +!! print *,'wake_dadv, k, xt1inv(1,k), xt1(1,k), xt1inv(1,k)-xt1(1,k) ', & +!! k, xt1inv(1,k), xt1(1,k), xt1inv(1,k)-xt1(1,k) +!! ENDIF +!! IF (abs(xt2inv(1,k)-xt2(1,k)) .GT. 1.e-15*xt2(1,k) ) THEN +!! print *,'wake_dadv, k, xt2inv(1,k), xt2(1,k), xt2inv(1,k)-xt2(1,k) ', & +!! k, xt2inv(1,k), xt2(1,k), xt2inv(1,k)-xt2(1,k) +!! ENDIF +!! IF (abs(xq1inv(1,k)-xq1(1,k)) .GT. 1.e-15*xq1(1,k) ) THEN +!! print *,'wake_dadv, k, xq1inv(1,k), xq1(1,k), xq1inv(1,k)-xq1(1,k) ', & +!! k, xq1inv(1,k), xq1(1,k), xq1inv(1,k)-xq1(1,k) +!! ENDIF +!! IF (abs(xq2inv(1,k)-xq2(1,k)) .GT. 1.e-15*xq2(1,k) ) THEN +!! print *,'wake_dadv, k, xq2inv(1,k), xq2(1,k), xq2inv(1,k)-xq2(1,k) ', & +!! k, xq2inv(1,k), xq2(1,k), xq2inv(1,k)-xq2(1,k) +!! ENDIF +!! ENDDO + + ELSE ! (flag_dadv_implicit) + +! Explicit scheme : compute directly d_deltat_dadv and d_tb_dadv +! (and similarly for d_deltaq_dadv and d_qb_dadv). +! ----------------------------------------------------------------------------------------- + + DO i = 1, klon + IF (wk_adv(i)) THEN !!! nrlmd + d_thx(i, 1) = 0. + d_thw(i, 1) = 0. + d_dth(i, 1) = 0. + d_qx(i, 1) = 0. + d_qw(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_thx(i, k) = thx(i, k-1) - thx(i, k) + d_thw(i, k) = thw(i, k-1) - thw(i, k) + d_dth(i, k) = delta_th(i, k-1) - delta_th(i, k) + d_qx(i, k) = qx(i, k-1) - qx(i, k) + d_qw(i, k) = qw(i, k-1) - qw(i, k) + d_dq(i, k) = delta_q(i, k-1) - delta_q(i, k) + END IF ! (wk_adv(i) .AND. k<=kupper(i)+1) + END DO + END DO + + DO k = 1, klev-1 + DO i = 1, klon + IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN + d_deltat_dadv(i, k) = dtime/(ph(i,k)-ph(i,k+1))* & + (rr22(i)*deltomg(i,k)*sigmaw(i)*d_thx(i,k) - & + rr21(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_thw(i,k+1) )*ppi(i, k) - & + dtime*entr_s(i,k)*delta_th(i,k)/sigmaw(i)*ppi(i, k) +! + d_deltaq_dadv(i, k) = dtime/(ph(i,k)-ph(i,k+1))* & + (rr22(i)*deltomg(i,k)*sigmaw(i)*d_qx(i,k)- & + rr21(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_qw(i,k+1) ) - & + dtime*entr_s(i,k)*delta_q(i,k)/sigmaw(i) + + ! and increment large scale tendencies + d_tb_dadv(i, k) = dtime*((rr12(i)*deltomg(i,k)*sigmaw(i)*d_thx(i,k)- & + rr11(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_thw(i,k+1))/ & + (ph(i,k)-ph(i,k+1)) & + -sigmaw(i)*(1.-sigmaw(i))*delta_th(i,k)*(deltomg(i,k)-deltomg(i,k+1))/ & + (ph(i,k)-ph(i,k+1)) )*ppi(i, k) + + d_qb_dadv(i, k) = dtime*((rr12(i)*deltomg(i,k)*sigmaw(i)*d_qx(i,k)- & + rr11(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_qw(i,k+1))/ & + (ph(i,k)-ph(i,k+1)) & + -sigmaw(i)*(1.-sigmaw(i))*delta_q(i,k)*(deltomg(i,k)-deltomg(i,k+1))/ & + (ph(i,k)-ph(i,k+1)) ) + ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN + d_tb_dadv(i, k) = dtime*(rr12(i)*deltomg(i, k)*sigmaw(i)*d_thx(i, k)/(ph(i, k)-ph(i, k+1)))*ppi(i, k) + + d_qb_dadv(i, k) = dtime*(rr12(i)*deltomg(i, k)*sigmaw(i)*d_qx(i, k)/(ph(i, k)-ph(i, k+1))) + END IF ! (wk_adv(i) .AND. k<=kupper(i)-1) + END DO + END DO + + ENDIF! (flag_dadv_implicit) + + END SUBROUTINE wake_dadv END MODULE lmdz_wake Index: LMDZ6/trunk/libf/phylmd/lmdz_wake_ini.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_wake_ini.f90 (revision 5760) +++ LMDZ6/trunk/libf/phylmd/lmdz_wake_ini.f90 (revision 5761) @@ -55,4 +55,7 @@ INTEGER, SAVE, PROTECTED :: iflag_wk_profile !$OMP THREADPRIVATE(iflag_wk_profile) + + LOGICAL, SAVE, PROTECTED :: flag_dadv_implicit + !$OMP THREADPRIVATE(flag_dadv_implicit) INTEGER, SAVE, PROTECTED :: wk_nsub @@ -250,4 +253,7 @@ CALL getin_p('wk_frac_int_delta_t', wk_frac_int_delta_t) + flag_dadv_implicit = .FALSE. + CALL getin_p('flag_dadv_implicit', flag_dadv_implicit) + CALL getin_p('CPPKEY_IOPHYS_WK', CPPKEY_IOPHYS_WK)