MODULE climb_hq_mod ! ! Module to solve the verctical diffusion of "q" and "H"; ! specific humidity and potential energi. ! USE dimphy #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso ! ajout C Risi pour isos #endif IMPLICIT NONE PRIVATE PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah !$OMP THREADPRIVATE(gamaq,gamah) REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q !$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q) REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H !$OMP THREADPRIVATE(Ccoef_H, Dcoef_H) REAL, DIMENSION(:), ALLOCATABLE :: Acoef_Q, Bcoef_Q !$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q) REAL, DIMENSION(:), ALLOCATABLE :: Acoef_H, Bcoef_H !$OMP THREADPRIVATE(Acoef_H, Bcoef_H) REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq !$OMP THREADPRIVATE(Kcoefhq) REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: h_old ! for diagnostics, h before solving diffusion !$OMP THREADPRIVATE(h_old) REAL, SAVE, DIMENSION(:), ALLOCATABLE :: d_h_col_vdf ! for diagnostics, vertical integral of enthalpy change !$OMP THREADPRIVATE(d_h_col_vdf) REAL, SAVE, DIMENSION(:), ALLOCATABLE :: f_h_bnd ! for diagnostics, enthalpy flux at surface !$OMP THREADPRIVATE(f_h_bnd) #ifdef ISO REAL, DIMENSION(:,:,:), ALLOCATABLE :: gamaxt !$OMP THREADPRIVATE(gamaxt) REAL, DIMENSION(:,:,:), ALLOCATABLE :: Ccoef_XT, Dcoef_XT !$OMP THREADPRIVATE(Ccoef_XT, Dcoef_XT) REAL, DIMENSION(:,:), ALLOCATABLE :: Acoef_XT, Bcoef_XT !$OMP THREADPRIVATE(Acoef_XT, Bcoef_XT) #endif CONTAINS ! !**************************************************************************************** ! SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, & delp, temp, q, dtime, & !!! nrlmd le 02/05/2011 Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, & Kcoef_hq_out, gama_q_out, gama_h_out, & !!! Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out & #ifdef ISO ,xt, & Ccoef_XT_out,Dcoef_XT_out,gama_xt_out, & Acoef_XT_out, Bcoef_XT_out & #endif ) #ifdef ISOVERIF USE isotopes_mod, ONLY: iso_eau,iso_HDO !USE isotopes_verif_mod, ONLY: errmax, errmaxrel USE isotopes_verif_mod #endif USE yomcst_mod_h USE compbl_mod_h ! This routine calculates recursivly the coefficients C and D ! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is ! the index of the vertical layer. ! ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: knon REAL, DIMENSION(klon,klev), INTENT(IN) :: coefhq REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs REAL, DIMENSION(klon,klev), INTENT(IN) :: temp, delp ! temperature REAL, DIMENSION(klon,klev), INTENT(IN) :: q REAL, INTENT(IN) :: dtime #ifdef ISO REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: xt #endif ! Output arguments !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_H_out REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_Q_out REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_H_out REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_Q_out #ifdef ISO REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: Acoef_XT_out REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: Bcoef_XT_out #endif !!! nrlmd le 02/05/2011 REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_H_out REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_Q_out REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_H_out REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_Q_out REAL, DIMENSION(klon,klev), INTENT(OUT) :: Kcoef_hq_out REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_q_out REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_h_out #ifdef ISO REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: Ccoef_XT_out REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: Dcoef_XT_out REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: gama_xt_out #endif !!! ! Local variables !**************************************************************************************** LOGICAL, SAVE :: first=.TRUE. !$OMP THREADPRIVATE(first) ! JLD now renamed h_old and declared in module ! REAL, DIMENSION(klon,klev) :: local_H REAL, DIMENSION(klon) :: psref REAL :: delz, pkh INTEGER :: k, i, ierr #ifdef ISO real, DIMENSION(klon,2:klev) :: gamaxt_tmp real, DIMENSION(klon,klev) :: xt_tmp real, DIMENSION(klon,klev) :: Ccoef_XT_tmp,Dcoef_XT_tmp real, DIMENSION(klon) :: Acoef_XT_tmp,Bcoef_XT_tmp integer ixt #endif #ifdef ISO #ifdef ISOVERIF if (iso_eau.gt.0) then do k = 1, klev DO i = 1, knon call iso_verif_egalite_choix( & xt(iso_eau,i,k),q(i,k), & 'climb_hq 100',errmax,errmaxrel) enddo enddo endif ! if (iso_eau.gt.0) then #endif #endif !**************************************************************************************** ! 1) ! Allocation at first time step only ! !**************************************************************************************** IF (first) THEN first=.FALSE. ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr ALLOCATE(Ccoef_H(klon,klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_H, ierr=', ierr ALLOCATE(Dcoef_H(klon,klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_H, ierr=', ierr ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr ALLOCATE(Kcoefhq(klon,klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc Kcoefhq, ierr=', ierr ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr ALLOCATE(gamah(1:klon,2:klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr #ifdef ISO ALLOCATE(Ccoef_XT(ntraciso,klon,klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_XT, ierr=', ierr ALLOCATE(Dcoef_XT(ntraciso,klon,klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_XT, ierr=', ierr ALLOCATE(Acoef_XT(ntraciso,klon), Bcoef_XT(ntraciso,klon), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_XT and Bcoef_XT, ierr=', ierr ALLOCATE(gamaxt(ntraciso,1:klon,2:klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaxt, ierr=', ierr #endif ALLOCATE(h_old(klon,klev), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc h_old, ierr=', ierr ALLOCATE(d_h_col_vdf(klon), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr ALLOCATE(f_h_bnd(klon), STAT=ierr) IF ( ierr /= 0 ) PRINT*,' pb in allloc f_h_bnd, ierr=', ierr END IF !**************************************************************************************** ! 2) ! Definition of the coeficient K ! !**************************************************************************************** Kcoefhq(:,:) = 0.0 DO k = 2, klev DO i = 1, knon Kcoefhq(i,k) = & coefhq(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) & *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2 ENDDO ENDDO !**************************************************************************************** ! 3) ! Calculation of gama for "Q" and "H" ! !**************************************************************************************** ! surface pressure is used as reference psref(:) = paprs(:,1) ! definition of gama IF (iflag_pbl == 1) THEN gamaq(:,:) = 0.0 gamah(:,:) = -1.0e-03 gamah(:,2) = -2.5e-03 #ifdef ISO do ixt=1,ntraciso gamaxt(:,:,:) = 0.0 enddo ! do ixt=1,ntraciso #endif ! conversion de gama DO k = 2, klev DO i = 1, knon delz = RD * (temp(i,k-1)+temp(i,k)) / & 2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k)) pkh = (psref(i)/paprs(i,k))**RKAPPA ! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches gamaq(i,k) = gamaq(i,k) * delz ! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches gamah(i,k) = gamah(i,k) * delz * RCPD * pkh #ifdef ISO do ixt=1,ntraciso gamaxt(ixt,i,k) = gamaxt(ixt,i,k) * delz enddo #endif ENDDO ENDDO ELSE gamaq(:,:) = 0.0 gamah(:,:) = 0.0 #ifdef ISO do ixt=1,ntraciso gamaxt(:,:,:) = 0.0 enddo ! do ixt=1,ntraciso #endif ENDIF #ifdef ISO #ifdef ISOVERIF do k = 2, klev DO i = 1, knon call iso_verif_egalite_choix( & gamaxt(iso_eau,i,k),gamaq(i,k), & 'climb_hq 209',errmax,errmaxrel) enddo enddo #endif #endif !**************************************************************************************** ! 4) ! Calculte the coefficients C and D for specific humidity, q ! !**************************************************************************************** CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), & Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q) #ifdef ISO do ixt=1,ntraciso ! compression do k = 2, klev DO i = 1, knon gamaxt_tmp(i,k)=gamaxt(ixt,i,k) enddo enddo !do k = 2, klev do k = 1, klev DO i = 1, knon xt_tmp(i,k)=xt(ixt,i,k) enddo enddo !do k = 2, klev !appel routine generique CALL calc_coef(knon, Kcoefhq(:,:), gamaxt_tmp(:,:), delp(:,:), xt_tmp(:,:), & Ccoef_XT_tmp(:,:), Dcoef_XT_tmp(:,:), Acoef_XT_tmp, Bcoef_XT_tmp) ! decompression do k = 1, klev DO i = 1, knon Ccoef_XT(ixt,i,k)=Ccoef_XT_tmp(i,k) Dcoef_XT(ixt,i,k)=Dcoef_XT_tmp(i,k) enddo enddo !do k = 2, klev DO i = 1, knon Acoef_XT(ixt,i)=Acoef_XT_tmp(i) Bcoef_XT(ixt,i)=Bcoef_XT_tmp(i) enddo enddo ! do ixt=1,ntraciso #ifdef ISOVERIF if (iso_eau.gt.0) then do k = 1, klev DO i = 1, knon call iso_verif_egalite_choix( & Ccoef_XT(iso_eau,i,k),Ccoef_Q(i,k), & 'climb_hq 234c',errmax,errmaxrel) call iso_verif_egalite_choix( & Dcoef_XT(iso_eau,i,k),Dcoef_Q(i,k), & 'climb_hq 234d',errmax,errmaxrel) enddo !DO i = 1, knon enddo !do k = 2, klev DO i = 1, knon call iso_verif_egalite_choix( & Acoef_XT(iso_eau,i),Acoef_Q(i), & 'climb_hq 234a',errmax,errmaxrel) call iso_verif_egalite_choix( & Bcoef_XT(iso_eau,i),Bcoef_Q(i), & 'climb_hq 234b',errmax,errmaxrel) enddo !DO i = 1, knon endif !if (iso_eau.gt.0) then #endif #endif !**************************************************************************************** ! 5) ! Calculte the coefficients C and D for potentiel entalpie, H ! !**************************************************************************************** h_old(:,:) = 0.0 DO k=1,klev DO i = 1, knon ! convertie la temperature en entalpie potentielle h_old(i,k) = RCPD * temp(i,k) * & (psref(i)/pplay(i,k))**RKAPPA ENDDO ENDDO CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), h_old(:,:), & Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H) !**************************************************************************************** ! 6) ! Return the first layer in output variables ! !**************************************************************************************** Acoef_H_out = Acoef_H Bcoef_H_out = Bcoef_H Acoef_Q_out = Acoef_Q Bcoef_Q_out = Bcoef_Q #ifdef ISO Acoef_XT_out = Acoef_XT Bcoef_XT_out = Bcoef_XT #endif !**************************************************************************************** ! 7) ! If Pbl is split, return also the other layers in output variables ! !**************************************************************************************** !!! jyg le 07/02/2012 !!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN IF (mod(iflag_pbl_split,10) .ge.1) THEN !!! nrlmd le 02/05/2011 DO k= 1, klev DO i= 1, klon Ccoef_H_out(i,k) = Ccoef_H(i,k) Dcoef_H_out(i,k) = Dcoef_H(i,k) Ccoef_Q_out(i,k) = Ccoef_Q(i,k) Dcoef_Q_out(i,k) = Dcoef_Q(i,k) Kcoef_hq_out(i,k) = Kcoefhq(i,k) #ifdef ISO do ixt=1,ntraciso Ccoef_XT_out(ixt,i,k) = Ccoef_XT(ixt,i,k) Dcoef_XT_out(ixt,i,k) = Dcoef_XT(ixt,i,k) enddo #endif IF (k.eq.1) THEN gama_h_out(i,k) = 0. gama_q_out(i,k) = 0. #ifdef ISO do ixt=1,ntraciso gama_xt_out(ixt,i,k) = 0. enddo #endif ELSE gama_h_out(i,k) = gamah(i,k) gama_q_out(i,k) = gamaq(i,k) #ifdef ISO do ixt=1,ntraciso gama_xt_out(ixt,i,k) = gamaxt(ixt,i,k) enddo #endif ENDIF ENDDO ENDDO !!! ENDIF ! (mod(iflag_pbl_split,2) .ge.1) !!! END SUBROUTINE climb_hq_down ! !**************************************************************************************** ! SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef) ! ! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1) ! where X is H or Q, and k the vertical level k=1,klev ! USE yomcst_mod_h ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: knon REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, delp REAL, DIMENSION(klon,klev), INTENT(IN) :: X REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama ! Output arguments !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef ! Local variables !**************************************************************************************** INTEGER :: k, i REAL :: buf !**************************************************************************************** ! Niveau au sommet, k=klev ! !**************************************************************************************** Ccoef(:,:) = 0.0 Dcoef(:,:) = 0.0 DO i = 1, knon buf = delp(i,klev) + Kcoef(i,klev) Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf Dcoef(i,klev) = Kcoef(i,klev)/buf END DO !**************************************************************************************** ! Niveau (klev-1) <= k <= 2 ! !**************************************************************************************** DO k=(klev-1),2,-1 DO i = 1, knon buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1)) Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + & Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf Dcoef(i,k) = Kcoef(i,k)/buf END DO END DO !**************************************************************************************** ! Niveau k=1 ! !**************************************************************************************** DO i = 1, knon buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2)) Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf Bcoef(i) = -1. * RG / buf END DO END SUBROUTINE calc_coef ! !**************************************************************************************** ! SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, & flx_q1, flx_h1, paprs, pplay, & !!! nrlmd le 02/05/2011 Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, & Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, & Kcoef_hq_in, gama_q_in, gama_h_in, & !!! flux_q, flux_h, d_q, d_t & #ifdef ISO ,xt_old, flx_xt1, & Acoef_XT_in,Bcoef_XT_in,Ccoef_XT_in,Dcoef_XT_in,gama_xt_in, & flux_xt, d_xt & #endif ) #ifdef ISOVERIF USE infotrac_phy, ONLY: nzone USE isotopes_mod, ONLY: iso_eau,iso_HDO,iso_O18, ridicule USE isotopes_verif_mod #endif USE yomcst_mod_h USE compbl_mod_h ! ! This routine calculates the flux and tendency of the specific humidity q and ! the potential engergi H. ! The quantities q and H are calculated according to ! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients ! C and D are known from before and k is index of the vertical layer. ! ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: knon REAL, INTENT(IN) :: dtime REAL, DIMENSION(klon,klev), INTENT(IN) :: t_old, q_old REAL, DIMENSION(klon), INTENT(IN) :: flx_q1, flx_h1 REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay !!! nrlmd le 02/05/2011 REAL, DIMENSION(klon), INTENT(IN) :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in REAL, DIMENSION(klon,klev), INTENT(IN) :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef_hq_in, gama_q_in, gama_h_in #ifdef ISO REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: Acoef_XT_in, Bcoef_XT_in REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: Ccoef_XT_in, Dcoef_XT_in REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: gama_xt_in #endif !!! ! Output arguments !**************************************************************************************** REAL, DIMENSION(klon,klev), INTENT(OUT) :: flux_q, flux_h, d_q, d_t ! Local variables !**************************************************************************************** LOGICAL, SAVE :: last=.FALSE. !$OMP THREADPRIVATE(last) REAL, DIMENSION(klon,klev) :: h_new, q_new REAL, DIMENSION(klon) :: psref INTEGER :: k, i, ierr #ifdef ISO REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: xt_old REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: flux_xt, d_xt REAL, DIMENSION(ntraciso,klon,klev) :: xt_new REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: flx_xt1 integer ixt !#ifdef ISOVERIF ! integer iso_verif_noNaN_nostop !#endif #endif !**************************************************************************************** ! 1) ! Definition of some variables REAL, DIMENSION(klon,klev) :: d_h, zairm ! !**************************************************************************************** #ifdef ISO #ifdef ISOVERIF DO k = 1, klev DO i = 1, knon if (iso_eau.gt.0) then call iso_verif_egalite(xt_old(iso_eau,i,k), & & q_old(i,k),'climb_hq_mod 421') endif #ifdef ISOTRAC IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 422') #endif enddo enddo #endif #endif flux_q(:,:) = 0.0 flux_h(:,:) = 0.0 d_q(:,:) = 0.0 d_t(:,:) = 0.0 d_h(:,:) = 0.0 f_h_bnd(:)= 0.0 #ifdef ISO q_new(:,:) = 0.0 flux_xt(:,:,:) = 0.0 d_xt(:,:,:) = 0.0 xt_new(:,:,:) = 0.0 #endif psref(1:knon) = paprs(1:knon,1) !!! jyg le 07/02/2012 !!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN IF (mod(iflag_pbl_split,10) .ge.1) THEN !!! nrlmd le 02/05/2011 DO i = 1, knon Acoef_H(i)=Acoef_H_in(i) Acoef_Q(i)=Acoef_Q_in(i) Bcoef_H(i)=Bcoef_H_in(i) Bcoef_Q(i)=Bcoef_Q_in(i) #ifdef ISO do ixt=1,ntraciso Acoef_XT(ixt,i)=Acoef_XT_in(ixt,i) Bcoef_XT(ixt,i)=Bcoef_XT_in(ixt,i) enddo #endif ENDDO DO k = 1, klev DO i = 1, knon Ccoef_H(i,k)=Ccoef_H_in(i,k) Ccoef_Q(i,k)=Ccoef_Q_in(i,k) Dcoef_H(i,k)=Dcoef_H_in(i,k) Dcoef_Q(i,k)=Dcoef_Q_in(i,k) Kcoefhq(i,k)=Kcoef_hq_in(i,k) #ifdef ISO do ixt=1,ntraciso Ccoef_XT(ixt,i,k)=Ccoef_XT_in(ixt,i,k) Dcoef_XT(ixt,i,k)=Dcoef_XT_in(ixt,i,k) enddo #endif IF (k.gt.1) THEN gamah(i,k)=gama_h_in(i,k) gamaq(i,k)=gama_q_in(i,k) #ifdef ISO do ixt=1,ntraciso gamaxt(ixt,i,k)=gama_xt_in(ixt,i,k) enddo #endif ENDIF ENDDO ENDDO !!! ENDIF ! (mod(iflag_pbl_split,2) .ge.1) !!! !**************************************************************************************** ! 2) ! Calculation of Q and H ! !**************************************************************************************** !- First layer q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)*flx_q1(1:knon)*dtime h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)*flx_h1(1:knon)*dtime f_h_bnd(1:knon) = flx_h1(1:knon) #ifdef ISO do ixt=1,ntraciso xt_new(ixt,1:knon,1) = Acoef_XT(ixt,1:knon) + Bcoef_XT(ixt,1:knon)*flx_xt1(ixt,1:knon)*dtime enddo ! do ixt=1,ntraciso #endif !- All the other layers DO k = 2, klev DO i = 1, knon q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1) h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1) #ifdef ISO do ixt=1,ntraciso xt_new(ixt,i,k) = Ccoef_XT(ixt,i,k) + Dcoef_XT(ixt,i,k)*xt_new(ixt,i,k-1) enddo ! do ixt=1,ntraciso #endif END DO END DO #ifdef ISO #ifdef ISOVERIF DO k = 1, klev DO i = 1, knon do ixt=1,ntraciso if (iso_verif_noNaN_nostop(xt_new(ixt,i,k),'climb_hq 507').eq.1) then write(*,*) 'Acoef_XT(ixt,i)=',Acoef_XT(ixt,i) write(*,*) 'Bcoef_XT(ixt,i)=',Bcoef_XT(ixt,i) write(*,*) 'flx_xt1(ixt,i)=',flx_xt1(ixt,i) if (k.ge.2) then write(*,*) 'Ccoef_XT(ixt,i,k)=',Ccoef_XT(ixt,i,k) write(*,*) 'Dcoef_XT(ixt,i,k)=',Dcoef_XT(ixt,i,k) endif stop endif enddo !do ixt=1,ntraciso END DO END DO #endif #ifdef ISOVERIF if (iso_eau.gt.0) then call iso_verif_egalite_vect2D( & & xt_new,q_new, & & 'climb_hq_mod 504',ntraciso,klon,klev) endif !if (iso_eau.gt.0) then if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then do i=1,klon do k=1,klev if (q_new(i,k).gt.ridicule) then if (iso_verif_o18_aberrant_nostop( & & xt_new(iso_HDO,i,k)/q_new(i,k), & & xt_new(iso_O18,i,k)/q_new(i,k), & & 'climb_hq_mod 690').eq.1) then write(*,*) 'i,k,q_new(i,k)=',i,k,q_new(i,k) stop endif ! if (iso_verif_o18_aberrant_nostop endif !if (q_seri(i,k).gt.errmax) then enddo !k=1,klev enddo !i=1,klon endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then #endif #endif !**************************************************************************************** ! 3) ! Calculation of the flux for Q and H ! !**************************************************************************************** !- The flux at first layer, k=1 flux_q(1:knon,1)=flx_q1(1:knon) flux_h(1:knon,1)=flx_h1(1:knon) #ifdef ISO do ixt=1,ntraciso flux_xt(ixt,1:knon,1)=flx_xt1(ixt,1:knon) enddo ! do ixt=1,ntraciso #endif !- The flux at all layers above surface DO k = 2, klev DO i = 1, knon flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * & (q_new(i,k)-q_new(i,k-1)+gamaq(i,k)) flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * & (h_new(i,k)-h_new(i,k-1)+gamah(i,k)) #ifdef ISO do ixt=1,ntraciso flux_xt(ixt,i,k) = (Kcoefhq(i,k)/RG/dtime) * & (xt_new(ixt,i,k)-xt_new(ixt,i,k-1)+gamaxt(ixt,i,k)) enddo ! do ixt=1,ntraciso #endif END DO END DO !**************************************************************************************** ! 4) ! Calculation of tendency for Q and H ! !**************************************************************************************** d_h_col_vdf(:) = 0.0 DO k = 1, klev DO i = 1, knon d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k) d_q(i,k) = q_new(i,k) - q_old(i,k) d_h(i,k) = h_new(i,k) - h_old(i,k) !JLD d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD !correction a venir ! layer air mass zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg d_h_col_vdf(i) = d_h_col_vdf(i) + d_h(i,k)*zairm(i,k) #ifdef ISO do ixt=1,ntraciso d_xt(ixt,i,k) = xt_new(ixt,i,k) - xt_old(ixt,i,k) enddo ! do ixt=1,ntraciso #ifdef ISOVERIF do ixt=1,ntraciso call iso_verif_noNaN(xt_new(ixt,i,k),'climb_hq 562') call iso_verif_noNaN(d_xt(ixt,i,k),'climb_hq 562') enddo ! do ixt=1,ntraciso #endif #ifdef ISOVERIF if (iso_eau.gt.0) then call iso_verif_egalite(d_xt(iso_eau,i,k), & & d_q(i,k),'climb_hq_mod 503') call iso_verif_egalite(xt_new(iso_eau,i,k), & & q_new(i,k),'climb_hq_mod 503b') endif #ifdef ISOTRAC IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 526') #endif #endif #endif END DO END DO #ifdef ISO #ifdef ISOVERIF ! write(*,*) 'climb_hq_mod 758: d_xt,d_q=',d_xt(iso_eau,1,1),d_q(1,1) if (iso_eau.gt.0) then call iso_verif_egalite_vect2D( & d_xt,d_q, & 'climb_hq_mod 761',ntraciso,klon,klev) endif #endif #endif !**************************************************************************************** ! Some deallocations ! !**************************************************************************************** IF (last) THEN DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr) IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr) IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr DEALLOCATE(gamaq, gamah,stat=ierr) IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr DEALLOCATE(Kcoefhq,stat=ierr) IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr) IF ( ierr /= 0 ) PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr #ifdef ISO DEALLOCATE(Ccoef_XT, Dcoef_XT ,stat=ierr) IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_XT, Dcoef_XT, ierr=', ierr DEALLOCATE(Acoef_XT, Bcoef_XT, stat=ierr) IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_XT, Bcoef_XT, ierr=', ierr DEALLOCATE(gamaxt, stat=ierr) IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaxt, ierr=', ierr #endif END IF END SUBROUTINE climb_hq_up ! !**************************************************************************************** ! END MODULE climb_hq_mod