Ignore:
Timestamp:
Apr 13, 2017, 7:07:54 PM (8 years ago)
Author:
jyg
Message:

Changes in module add_phys_tend_mod.F90: (i)
vertical integrations are performed in subroutine
integr_v; (ii) the new subroutine diag_phys_tend
makes it possible to check water and energy
conservation of a set of tendencies added to any
atmospheric state.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/add_phys_tend_mod.F90

    r2812 r2848  
    223223
    224224  if (fl_ebil .GT. 0) then
    225   ! Reset variables
    226     zqw_col(:,:) = 0.
    227     zql_col(:,:) = 0.
    228     zqs_col(:,:) = 0.
    229     zek_col(:,:) = 0.
    230     zh_dair_col(:,:) = 0.
    231     zh_qw_col(:,:) = 0.
    232     zh_ql_col(:,:) = 0.
    233     zh_qs_col(:,:) = 0.
     225    ! ------------------------------------------------
     226    ! Compute vertical sum for each atmospheric column
     227    ! ------------------------------------------------
     228    n=1   ! begining of time step
    234229
    235230    zcpvap = rcpv
    236231    zcwat = rcw
    237232    zcice = rcs
    238 !JLD    write (*,*) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice
    239  
    240     ! ------------------------------------------------
    241     ! Compute vertical sum for each atmospheric column
    242     ! ------------------------------------------------
    243     n=1   ! begining of time step
    244     DO k = 1, klev
    245       DO i = 1, klon
    246         ! Watter mass
    247         zqw_col(i,n) = zqw_col(i,n) + q_seri(i, k)*zairm(i, k)
    248         zql_col(i,n) = zql_col(i,n) + ql_seri(i, k)*zairm(i, k)
    249         zqs_col(i,n) = zqs_col(i,n) + qs_seri(i, k)*zairm(i, k)
    250         ! Kinetic Energy
    251         zek_col(i,n) = zek_col(i,n) + 0.5*(u_seri(i,k)**2+v_seri(i,k)**2)*zairm(i, k)
    252         ! Air enthalpy : dry air, water vapour, liquid, solid
    253         zh_dair_col(i,n) = zh_dair_col(i,n) + rcpd*(1.-q_seri(i,k)-ql_seri(i,k)-qs_seri(i,k))* &
    254           zairm(i, k)*t_seri(i, k)
    255         zh_qw_col(i,n) = zh_qw_col(i,n) +  zcpvap*t_seri(i, k)         *q_seri(i, k)*zairm(i, k)    !jyg
    256         zh_ql_col(i,n) = zh_ql_col(i,n) + (zcpvap*t_seri(i, k) - rlvtt)*ql_seri(i, k)*zairm(i, k)   !jyg
    257         zh_qs_col(i,n) = zh_qs_col(i,n) + (zcpvap*t_seri(i, k) - rlstt)*qs_seri(i, k)*zairm(i, k)   !jyg
    258       END DO
    259     END DO
    260     ! compute total air enthalpy
    261     zh_col(:,n) = zh_dair_col(:,n) + zh_qw_col(:,n) + zh_ql_col(:,n) + zh_qs_col(:,n)
     233
     234    CALL integr_v(klon, klev, zcpvap, &
     235                  t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
     236                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
     237                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
    262238
    263239  end if ! end if (fl_ebil .GT. 0)
    … …  
    465441    ! ------------------------------------------------
    466442    n=2   ! end of time step
    467     DO k = 1, klev
    468       DO i = 1, klon
    469         ! Watter mass
    470         zqw_col(i,n) = zqw_col(i,n) + q_seri(i, k)*zairm(i, k)
    471         zql_col(i,n) = zql_col(i,n) + ql_seri(i, k)*zairm(i, k)
    472         zqs_col(i,n) = zqs_col(i,n) + qs_seri(i, k)*zairm(i, k)
    473         ! Kinetic Energy
    474         zek_col(i,n) = zek_col(i,n) + 0.5*(u_seri(i,k)**2+v_seri(i,k)**2)*zairm(i, k)
    475         ! Air enthalpy : dry air, water vapour, liquid, solid
    476         zh_dair_col(i,n) = zh_dair_col(i,n) + rcpd*(1.-q_seri(i,k)-ql_seri(i,k)-qs_seri(i,k))* &
    477           zairm(i, k)*t_seri(i, k)
    478         zh_qw_col(i,n) = zh_qw_col(i,n) +  zcpvap*t_seri(i, k)         *q_seri(i, k)*zairm(i, k)     !jyg
    479         zh_ql_col(i,n) = zh_ql_col(i,n) + (zcpvap*t_seri(i, k) - rlvtt)*ql_seri(i, k)*zairm(i, k)    !jyg
    480         zh_qs_col(i,n) = zh_qs_col(i,n) + (zcpvap*t_seri(i, k) - rlstt)*qs_seri(i, k)*zairm(i, k)    !jyg
    481       END DO
    482     END DO
    483     ! compute total air enthalpy
    484     zh_col(:,n) = zh_dair_col(:,n) + zh_qw_col(:,n) + zh_ql_col(:,n) + zh_qs_col(:,n)
     443
     444    CALL integr_v(klon, klev, zcpvap, &
     445                  t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
     446                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
     447                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
    485448
    486449    ! ------------------------------------------------
    … …  
    517480  RETURN
    518481END SUBROUTINE add_phys_tend
     482
     483SUBROUTINE diag_phys_tend (nlon, nlev, uu, vv, temp, qv, ql, qs, &
     484                          zdu,zdv,zdt,zdq,zdql,zdqs,paprs,text)
     485!======================================================================
     486! Ajoute les tendances des variables physiques aux variables
     487! d'etat de la dynamique t_seri, q_seri ...
     488! On en profite pour faire des tests sur les tendances en question.
     489!======================================================================
     490
     491
     492!======================================================================
     493! Declarations
     494!======================================================================
     495
     496USE phys_state_var_mod, ONLY : dtime, ftsol
     497USE geometry_mod, ONLY: longitude_deg, latitude_deg
     498USE print_control_mod, ONLY: prt_level
     499USE cmp_seri_mod
     500USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
     501  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
     502IMPLICIT none
     503  include "YOMCST.h"
     504  include "clesphys.h"
     505
     506! Arguments :
     507!------------
     508INTEGER, INTENT(IN)                           :: nlon, nlev
     509REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: uu, vv
     510REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: temp, qv, ql, qs
     511REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdu, zdv
     512REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdt, zdq, zdql, zdqs
     513REAL, DIMENSION(nlon,nlev+1),   INTENT(IN)    :: paprs
     514CHARACTER*(*),                  INTENT(IN)    :: text
     515
     516! Local :
     517!--------
     518REAL, DIMENSION(nlon,nlev)      :: uu_n, vv_n
     519REAL, DIMENSION(nlon,nlev)      :: temp_n, qv_n, ql_n, qs_n
     520
     521
     522!
     523INTEGER k, n
     524
     525integer debug_level
     526logical, save :: first=.true.
     527!$OMP THREADPRIVATE(first)
     528!
     529!======================================================================
     530! Variables for energy conservation tests
     531!======================================================================
     532!
     533
     534! zh_col-------  total enthalpy of vertical air column
     535! (air with watter vapour, liquid and solid) (J/m2)
     536! zh_dair_col--- total enthalpy of dry air (J/m2)
     537! zh_qw_col----  total enthalpy of watter vapour (J/m2)
     538! zh_ql_col----  total enthalpy of liquid watter (J/m2)
     539! zh_qs_col----  total enthalpy of solid watter  (J/m2)
     540! zqw_col------  total mass of watter vapour (kg/m2)
     541! zql_col------  total mass of liquid watter (kg/m2)
     542! zqs_col------  total mass of solid watter (kg/m2)
     543! zek_col------  total kinetic energy (kg/m2)
     544!
     545REAL zairm(nlon, nlev) ! layer air mass (kg/m2)
     546REAL zqw_col(nlon,2)
     547REAL zql_col(nlon,2)
     548REAL zqs_col(nlon,2)
     549REAL zek_col(nlon,2)
     550REAL zh_dair_col(nlon,2)
     551REAL zh_qw_col(nlon,2), zh_ql_col(nlon,2), zh_qs_col(nlon,2)
     552REAL zh_col(nlon,2)
     553
     554!======================================================================
     555! Initialisations
     556
     557     IF (prt_level >= 5) then
     558        write (*,*) "In diag_phys_tend, after ",text
     559        call flush
     560     end if
     561
     562     debug_level=10
     563     if (first) then
     564        print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs
     565        first=.false.
     566     endif
     567!
     568!  print *,'add_phys_tend: paprs ',paprs
     569!======================================================================
     570! Diagnostics for energy conservation tests
     571!======================================================================
     572  DO k = 1, nlev
     573    ! layer air mass
     574    zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg
     575  END DO
     576
     577  if (fl_ebil .GT. 0) then
     578    ! ------------------------------------------------
     579    ! Compute vertical sum for each atmospheric column
     580    ! ------------------------------------------------
     581    n=1   ! begining of time step
     582
     583    CALL integr_v(nlon, nlev, rcpv, &
     584                  temp, qv, ql, qs, uu, vv, zairm, &
     585                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
     586                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
     587
     588  end if ! end if (fl_ebil .GT. 0)
     589
     590!======================================================================
     591! Ajout des tendances sur le vent, la temperature et les diverses phases de l'eau
     592!======================================================================
     593
     594     uu_n(:,:)=uu(:,:)+zdu(:,:)
     595     vv_n(:,:)=vv(:,:)+zdv(:,:)
     596     qv_n(:,:)=qv(:,:)+zdq(:,:)
     597     ql_n(:,:)=ql(:,:)+zdql(:,:)
     598     qs_n(:,:)=qs(:,:)+zdqs(:,:)
     599     temp_n(:,:)=temp(:,:)+zdt(:,:)
     600
     601
     602
     603!======================================================================
     604! Diagnostics for energy conservation tests
     605!======================================================================
     606
     607  if (fl_ebil .GT. 0) then
     608 
     609    ! ------------------------------------------------
     610    ! Compute vertical sum for each atmospheric column
     611    ! ------------------------------------------------
     612    n=2   ! end of time step
     613
     614    CALL integr_v(nlon, nlev, rcpv, &
     615                  temp_n, qv_n, ql_n, qs_n, uu_n, vv_n, zairm, &
     616                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
     617                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
     618
     619    ! ------------------------------------------------
     620    ! Compute the changes by unit of time
     621    ! ------------------------------------------------
     622
     623    d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/dtime
     624    d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/dtime
     625    d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/dtime
     626    d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
     627
     628    d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/dtime
     629
     630   print *,'zdu ', zdu
     631   print *,'zdv ', zdv
     632   print *,'d_ek_col, zek_col(2), zek_col(1) ',d_ek_col(1), zek_col(1,2), zek_col(1,1)
     633
     634    d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/dtime
     635    d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/dtime
     636    d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/dtime
     637    d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/dtime
     638
     639    d_h_col = (zh_col(:,2)-zh_col(:,1))/dtime
     640
     641  end if ! end if (fl_ebil .GT. 0)
     642!
     643
     644  RETURN
     645END SUBROUTINE diag_phys_tend
     646
     647SUBROUTINE integr_v(nlon, nlev, zcpvap, &
     648                    temp, qv, ql, qs, uu, vv, zairm,  &
     649                    zqw_col, zql_col, zqs_col, zek_col, zh_dair_col, &
     650                    zh_qw_col, zh_ql_col, zh_qs_col, zh_col)
     651
     652IMPLICIT none
     653  include "YOMCST.h"
     654
     655INTEGER,                    INTENT(IN)    :: nlon,nlev
     656REAL,                       INTENT(IN)    :: zcpvap
     657REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: temp, qv, ql, qs, uu, vv
     658REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: zairm
     659REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqw_col
     660REAL, DIMENSION(nlon),      INTENT(OUT)   :: zql_col
     661REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqs_col
     662REAL, DIMENSION(nlon),      INTENT(OUT)   :: zek_col
     663REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_dair_col
     664REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qw_col
     665REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_ql_col
     666REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qs_col
     667REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_col
     668
     669INTEGER          :: i, k
     670
     671
     672  ! Reset variables
     673    zqw_col(:) = 0.
     674    zql_col(:) = 0.
     675    zqs_col(:) = 0.
     676    zek_col(:) = 0.
     677    zh_dair_col(:) = 0.
     678    zh_qw_col(:) = 0.
     679    zh_ql_col(:) = 0.
     680    zh_qs_col(:) = 0.
     681
     682!JLD    write (*,*) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice
     683 
     684    ! ------------------------------------------------
     685    ! Compute vertical sum for each atmospheric column
     686    ! ------------------------------------------------
     687    DO k = 1, nlev
     688      DO i = 1, nlon
     689        ! Watter mass
     690        zqw_col(i) = zqw_col(i) + qv(i, k)*zairm(i, k)
     691        zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
     692        zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
     693        ! Kinetic Energy
     694        zek_col(i) = zek_col(i) + 0.5*(uu(i,k)**2+vv(i,k)**2)*zairm(i, k)
     695        ! Air enthalpy : dry air, water vapour, liquid, solid
     696        zh_dair_col(i) = zh_dair_col(i) + rcpd*(1.-qv(i,k)-ql(i,k)-qs(i,k))* &
     697                                               zairm(i, k)*temp(i, k)
     698        zh_qw_col(i) = zh_qw_col(i) +  zcpvap*temp(i, k)         *qv(i, k)*zairm(i, k)    !jyg
     699        zh_ql_col(i) = zh_ql_col(i) + (zcpvap*temp(i, k) - rlvtt)*ql(i, k)*zairm(i, k)   !jyg
     700        zh_qs_col(i) = zh_qs_col(i) + (zcpvap*temp(i, k) - rlstt)*qs(i, k)*zairm(i, k)   !jyg
     701      END DO
     702    END DO
     703    ! compute total air enthalpy
     704    zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:)
     705
     706END SUBROUTINE integr_v
    519707
    520708SUBROUTINE prt_enerbil (text, itap)
    … …  
    624812      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)
    625813    end if
     814  CASE("convection") specific_diag
     815    if ( prt_level .GE. 5) then
     816      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)
     817      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)
     818    end if
    626819  END SELECT specific_diag
    627820
    628 9000 format (1x,A8,2x,A30,10E15.6)
     8219000 format (1x,A8,2x,A35,10E15.6)
    629822
    630823end if ! end if (fl_ebil .GT. 0)
Note: See TracChangeset for help on using the changeset viewer.