Index: /LMDZ5/trunk/libf/phylmd/add_phys_tend_mod.F90 =================================================================== --- /LMDZ5/trunk/libf/phylmd/add_phys_tend_mod.F90 (revision 2847) +++ /LMDZ5/trunk/libf/phylmd/add_phys_tend_mod.F90 (revision 2848) @@ -223,41 +223,17 @@ if (fl_ebil .GT. 0) then - ! Reset variables - zqw_col(:,:) = 0. - zql_col(:,:) = 0. - zqs_col(:,:) = 0. - zek_col(:,:) = 0. - zh_dair_col(:,:) = 0. - zh_qw_col(:,:) = 0. - zh_ql_col(:,:) = 0. - zh_qs_col(:,:) = 0. + ! ------------------------------------------------ + ! Compute vertical sum for each atmospheric column + ! ------------------------------------------------ + n=1 ! begining of time step zcpvap = rcpv zcwat = rcw zcice = rcs -!JLD write (*,*) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice - - ! ------------------------------------------------ - ! Compute vertical sum for each atmospheric column - ! ------------------------------------------------ - n=1 ! begining of time step - DO k = 1, klev - DO i = 1, klon - ! Watter mass - zqw_col(i,n) = zqw_col(i,n) + q_seri(i, k)*zairm(i, k) - zql_col(i,n) = zql_col(i,n) + ql_seri(i, k)*zairm(i, k) - zqs_col(i,n) = zqs_col(i,n) + qs_seri(i, k)*zairm(i, k) - ! Kinetic Energy - zek_col(i,n) = zek_col(i,n) + 0.5*(u_seri(i,k)**2+v_seri(i,k)**2)*zairm(i, k) - ! Air enthalpy : dry air, water vapour, liquid, solid - zh_dair_col(i,n) = zh_dair_col(i,n) + rcpd*(1.-q_seri(i,k)-ql_seri(i,k)-qs_seri(i,k))* & - zairm(i, k)*t_seri(i, k) - zh_qw_col(i,n) = zh_qw_col(i,n) + zcpvap*t_seri(i, k) *q_seri(i, k)*zairm(i, k) !jyg - zh_ql_col(i,n) = zh_ql_col(i,n) + (zcpvap*t_seri(i, k) - rlvtt)*ql_seri(i, k)*zairm(i, k) !jyg - zh_qs_col(i,n) = zh_qs_col(i,n) + (zcpvap*t_seri(i, k) - rlstt)*qs_seri(i, k)*zairm(i, k) !jyg - END DO - END DO - ! compute total air enthalpy - zh_col(:,n) = zh_dair_col(:,n) + zh_qw_col(:,n) + zh_ql_col(:,n) + zh_qs_col(:,n) + + CALL integr_v(klon, klev, zcpvap, & + t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, & + zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), & + zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n)) end if ! end if (fl_ebil .GT. 0) @@ -465,22 +441,9 @@ ! ------------------------------------------------ n=2 ! end of time step - DO k = 1, klev - DO i = 1, klon - ! Watter mass - zqw_col(i,n) = zqw_col(i,n) + q_seri(i, k)*zairm(i, k) - zql_col(i,n) = zql_col(i,n) + ql_seri(i, k)*zairm(i, k) - zqs_col(i,n) = zqs_col(i,n) + qs_seri(i, k)*zairm(i, k) - ! Kinetic Energy - zek_col(i,n) = zek_col(i,n) + 0.5*(u_seri(i,k)**2+v_seri(i,k)**2)*zairm(i, k) - ! Air enthalpy : dry air, water vapour, liquid, solid - zh_dair_col(i,n) = zh_dair_col(i,n) + rcpd*(1.-q_seri(i,k)-ql_seri(i,k)-qs_seri(i,k))* & - zairm(i, k)*t_seri(i, k) - zh_qw_col(i,n) = zh_qw_col(i,n) + zcpvap*t_seri(i, k) *q_seri(i, k)*zairm(i, k) !jyg - zh_ql_col(i,n) = zh_ql_col(i,n) + (zcpvap*t_seri(i, k) - rlvtt)*ql_seri(i, k)*zairm(i, k) !jyg - zh_qs_col(i,n) = zh_qs_col(i,n) + (zcpvap*t_seri(i, k) - rlstt)*qs_seri(i, k)*zairm(i, k) !jyg - END DO - END DO - ! compute total air enthalpy - zh_col(:,n) = zh_dair_col(:,n) + zh_qw_col(:,n) + zh_ql_col(:,n) + zh_qs_col(:,n) + + CALL integr_v(klon, klev, zcpvap, & + t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, & + zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), & + zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n)) ! ------------------------------------------------ @@ -517,4 +480,229 @@ RETURN END SUBROUTINE add_phys_tend + +SUBROUTINE diag_phys_tend (nlon, nlev, uu, vv, temp, qv, ql, qs, & + zdu,zdv,zdt,zdq,zdql,zdqs,paprs,text) +!====================================================================== +! Ajoute les tendances des variables physiques aux variables +! d'etat de la dynamique t_seri, q_seri ... +! On en profite pour faire des tests sur les tendances en question. +!====================================================================== + + +!====================================================================== +! Declarations +!====================================================================== + +USE phys_state_var_mod, ONLY : dtime, ftsol +USE geometry_mod, ONLY: longitude_deg, latitude_deg +USE print_control_mod, ONLY: prt_level +USE cmp_seri_mod +USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col & + & , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col +IMPLICIT none + include "YOMCST.h" + include "clesphys.h" + +! Arguments : +!------------ +INTEGER, INTENT(IN) :: nlon, nlev +REAL, DIMENSION(nlon,nlev), INTENT(IN) :: uu, vv +REAL, DIMENSION(nlon,nlev), INTENT(IN) :: temp, qv, ql, qs +REAL, DIMENSION(nlon,nlev), INTENT(IN) :: zdu, zdv +REAL, DIMENSION(nlon,nlev), INTENT(IN) :: zdt, zdq, zdql, zdqs +REAL, DIMENSION(nlon,nlev+1), INTENT(IN) :: paprs +CHARACTER*(*), INTENT(IN) :: text + +! Local : +!-------- +REAL, DIMENSION(nlon,nlev) :: uu_n, vv_n +REAL, DIMENSION(nlon,nlev) :: temp_n, qv_n, ql_n, qs_n + + +! +INTEGER k, n + +integer debug_level +logical, save :: first=.true. +!$OMP THREADPRIVATE(first) +! +!====================================================================== +! Variables for energy conservation tests +!====================================================================== +! + +! zh_col------- total enthalpy of vertical air column +! (air with watter vapour, liquid and solid) (J/m2) +! zh_dair_col--- total enthalpy of dry air (J/m2) +! zh_qw_col---- total enthalpy of watter vapour (J/m2) +! zh_ql_col---- total enthalpy of liquid watter (J/m2) +! zh_qs_col---- total enthalpy of solid watter (J/m2) +! zqw_col------ total mass of watter vapour (kg/m2) +! zql_col------ total mass of liquid watter (kg/m2) +! zqs_col------ total mass of solid watter (kg/m2) +! zek_col------ total kinetic energy (kg/m2) +! +REAL zairm(nlon, nlev) ! layer air mass (kg/m2) +REAL zqw_col(nlon,2) +REAL zql_col(nlon,2) +REAL zqs_col(nlon,2) +REAL zek_col(nlon,2) +REAL zh_dair_col(nlon,2) +REAL zh_qw_col(nlon,2), zh_ql_col(nlon,2), zh_qs_col(nlon,2) +REAL zh_col(nlon,2) + +!====================================================================== +! Initialisations + + IF (prt_level >= 5) then + write (*,*) "In diag_phys_tend, after ",text + call flush + end if + + debug_level=10 + if (first) then + print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs + first=.false. + endif +! +! print *,'add_phys_tend: paprs ',paprs +!====================================================================== +! Diagnostics for energy conservation tests +!====================================================================== + DO k = 1, nlev + ! layer air mass + zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg + END DO + + if (fl_ebil .GT. 0) then + ! ------------------------------------------------ + ! Compute vertical sum for each atmospheric column + ! ------------------------------------------------ + n=1 ! begining of time step + + CALL integr_v(nlon, nlev, rcpv, & + temp, qv, ql, qs, uu, vv, zairm, & + zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), & + zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n)) + + end if ! end if (fl_ebil .GT. 0) + +!====================================================================== +! Ajout des tendances sur le vent, la temperature et les diverses phases de l'eau +!====================================================================== + + uu_n(:,:)=uu(:,:)+zdu(:,:) + vv_n(:,:)=vv(:,:)+zdv(:,:) + qv_n(:,:)=qv(:,:)+zdq(:,:) + ql_n(:,:)=ql(:,:)+zdql(:,:) + qs_n(:,:)=qs(:,:)+zdqs(:,:) + temp_n(:,:)=temp(:,:)+zdt(:,:) + + + +!====================================================================== +! Diagnostics for energy conservation tests +!====================================================================== + + if (fl_ebil .GT. 0) then + + ! ------------------------------------------------ + ! Compute vertical sum for each atmospheric column + ! ------------------------------------------------ + n=2 ! end of time step + + CALL integr_v(nlon, nlev, rcpv, & + temp_n, qv_n, ql_n, qs_n, uu_n, vv_n, zairm, & + zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), & + zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n)) + + ! ------------------------------------------------ + ! Compute the changes by unit of time + ! ------------------------------------------------ + + d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/dtime + d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/dtime + d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/dtime + d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:) + + d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/dtime + + print *,'zdu ', zdu + print *,'zdv ', zdv + print *,'d_ek_col, zek_col(2), zek_col(1) ',d_ek_col(1), zek_col(1,2), zek_col(1,1) + + d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/dtime + d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/dtime + d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/dtime + d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/dtime + + d_h_col = (zh_col(:,2)-zh_col(:,1))/dtime + + end if ! end if (fl_ebil .GT. 0) +! + + RETURN +END SUBROUTINE diag_phys_tend + +SUBROUTINE integr_v(nlon, nlev, zcpvap, & + temp, qv, ql, qs, uu, vv, zairm, & + zqw_col, zql_col, zqs_col, zek_col, zh_dair_col, & + zh_qw_col, zh_ql_col, zh_qs_col, zh_col) + +IMPLICIT none + include "YOMCST.h" + +INTEGER, INTENT(IN) :: nlon,nlev +REAL, INTENT(IN) :: zcpvap +REAL, DIMENSION(nlon,nlev), INTENT(IN) :: temp, qv, ql, qs, uu, vv +REAL, DIMENSION(nlon,nlev), INTENT(IN) :: zairm +REAL, DIMENSION(nlon), INTENT(OUT) :: zqw_col +REAL, DIMENSION(nlon), INTENT(OUT) :: zql_col +REAL, DIMENSION(nlon), INTENT(OUT) :: zqs_col +REAL, DIMENSION(nlon), INTENT(OUT) :: zek_col +REAL, DIMENSION(nlon), INTENT(OUT) :: zh_dair_col +REAL, DIMENSION(nlon), INTENT(OUT) :: zh_qw_col +REAL, DIMENSION(nlon), INTENT(OUT) :: zh_ql_col +REAL, DIMENSION(nlon), INTENT(OUT) :: zh_qs_col +REAL, DIMENSION(nlon), INTENT(OUT) :: zh_col + +INTEGER :: i, k + + + ! Reset variables + zqw_col(:) = 0. + zql_col(:) = 0. + zqs_col(:) = 0. + zek_col(:) = 0. + zh_dair_col(:) = 0. + zh_qw_col(:) = 0. + zh_ql_col(:) = 0. + zh_qs_col(:) = 0. + +!JLD write (*,*) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice + + ! ------------------------------------------------ + ! Compute vertical sum for each atmospheric column + ! ------------------------------------------------ + DO k = 1, nlev + DO i = 1, nlon + ! Watter mass + zqw_col(i) = zqw_col(i) + qv(i, k)*zairm(i, k) + zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k) + zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k) + ! Kinetic Energy + zek_col(i) = zek_col(i) + 0.5*(uu(i,k)**2+vv(i,k)**2)*zairm(i, k) + ! Air enthalpy : dry air, water vapour, liquid, solid + zh_dair_col(i) = zh_dair_col(i) + rcpd*(1.-qv(i,k)-ql(i,k)-qs(i,k))* & + zairm(i, k)*temp(i, k) + zh_qw_col(i) = zh_qw_col(i) + zcpvap*temp(i, k) *qv(i, k)*zairm(i, k) !jyg + zh_ql_col(i) = zh_ql_col(i) + (zcpvap*temp(i, k) - rlvtt)*ql(i, k)*zairm(i, k) !jyg + zh_qs_col(i) = zh_qs_col(i) + (zcpvap*temp(i, k) - rlstt)*qs(i, k)*zairm(i, k) !jyg + END DO + END DO + ! compute total air enthalpy + zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:) + +END SUBROUTINE integr_v SUBROUTINE prt_enerbil (text, itap) @@ -624,7 +812,12 @@ write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_lsc(1), rlstt * snow_lsc(1), -(rcs-rcpd)*t_seri(1,1) * snow_lsc(1) end if + CASE("convection") specific_diag + if ( prt_level .GE. 5) then + write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_con(1), rlvtt * rain_con(1), -(rcw-rcpd)*t_seri(1,1) * rain_con(1) + write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_con(1), rlstt * snow_con(1), -(rcs-rcpd)*t_seri(1,1) * snow_con(1) + end if END SELECT specific_diag -9000 format (1x,A8,2x,A30,10E15.6) +9000 format (1x,A8,2x,A35,10E15.6) end if ! end if (fl_ebil .GT. 0)