LCOV - code coverage report
Current view: top level - dynamics/se - dycore_budget.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 4 329 1.2 %
Date: 2024-12-17 17:57:11 Functions: 1 2 50.0 %

          Line data    Source code
       1             : module dycore_budget
       2             : use shr_kind_mod, only: r8=>shr_kind_r8
       3             : implicit none
       4             : 
       5             : public :: print_budget
       6             : real(r8), parameter :: eps      = 1.0E-7_r8
       7             : real(r8), parameter :: eps_mass = 1.0E-12_r8
       8             : 
       9             : real(r8), save :: previous_dEdt_adiabatic_dycore    = 0.0_r8
      10             : real(r8), save :: previous_dEdt_dry_mass_adjust     = 0.0_r8
      11             : real(r8), save :: previous_dEdt_phys_dyn_coupl_err  = 0.0_r8
      12             : !=========================================================================================
      13             : contains
      14             : !=========================================================================================
      15             : 
      16      369408 : subroutine print_budget(hstwr)
      17             : 
      18             :   use spmd_utils,             only: masterproc
      19             :   use cam_abortutils,         only: endrun  
      20             :   use cam_logfile,            only: iulog
      21             :   use cam_budget,             only: cam_budget_get_global, is_cam_budget, thermo_budget_histfile_num, thermo_budget_history
      22             :   use cam_thermo,             only: thermo_budget_vars_descriptor, thermo_budget_num_vars, thermo_budget_vars_massv, &
      23             :                                     teidx, seidx, keidx, poidx
      24             :   use dimensions_mod,         only: use_cslam
      25             :   use control_mod,            only: ftype
      26             : 
      27             :   ! arguments
      28             :   logical, intent(in) :: hstwr(:)
      29             : 
      30             :   ! Local variables
      31             :   character(len=*), parameter :: subname = 'dycore_budget:print_budgets:'
      32             :   !
      33             :   ! physics energy tendencies
      34             :   !
      35             :   integer  :: idx(4)
      36             :   real(r8) :: dEdt_param_physE(4)      ! dE/dt CAM physics using physics E formula (phAP-phBP)
      37             :   real(r8) :: dEdt_param_dynE(4)       ! dE/dt CAM physics using dycore E (dyAP-dyBP)
      38             : 
      39             :   real(r8) :: dEdt_efix_physE(4)       ! dE/dt energy fixer using physics E formula (phBP-phBF)
      40             :   real(r8) :: dEdt_efix_dynE(4)        ! dE/dt energy fixer using dycore E formula (dyBP-dyBF)
      41             : 
      42             :   real(r8) :: dEdt_dme_adjust_physE(4) ! dE/dt dry mass adjustment using physics E formula (phAM-phAP)
      43             :   real(r8) :: dEdt_dme_adjust_dynE(4)  ! dE/dt dry mass adjustment using dycore E (dyAM-dyAP)
      44             : 
      45             :   real(r8) :: dEdt_param_efix_physE(4) ! dE/dt CAM physics + energy fixer using physics E formula (phAP-phBF)
      46             :   real(r8) :: dEdt_param_efix_dynE(4)  ! dE/dt CAM physics + energy fixer using dycore E formula (dyAP-dyBF)
      47             : 
      48             :   real(r8) :: dEdt_phys_total_dynE(4)  ! dE/dt physics total using dycore E (dyAM-dyBF)
      49             :                                        ! physics total = parameterizations + efix + dry-mass adjustment
      50             :   !
      51             :   ! SE dycore specific energy tendencies
      52             :   !
      53             :   real(r8) :: dEdt_phys_total_in_dyn(4) ! dEdt of physics total in dynamical core
      54             :                                         ! physics total = parameterizations + efix + dry-mass adjustment
      55             :   real(r8) :: dEdt_dycore_phys          ! dEdt dycore (estimated in physics)
      56             :   !
      57             :   ! mass budgets physics
      58             :   !
      59             :   real(r8) :: dMdt_efix                 ! mass tendency energy fixer
      60             :   real(r8) :: dMdt_parameterizations    ! mass tendency physics paramterizations
      61             :   real(r8) :: dMdt_dme_adjust           ! mass tendency dry-mass adjustment
      62             :   real(r8) :: dMdt_phys_total           ! mass tendency physics total (energy fixer + parameterizations + dry-mass adjustment)
      63             :   !
      64             :   ! mass budgets dynamics
      65             :   !
      66             :   real(r8) :: dMdt_floating_dyn         ! mass tendency floating dynamics (dAL-dBL)
      67             :   real(r8) :: dMdt_vert_remap           ! mass tendency vertical remapping (dAR-dAD)
      68             :   real(r8) :: dMdt_del4_fric_heat       ! mass tendency del4 frictional heating (dAH-dCH)
      69             :   real(r8) :: dMdt_del4_tot             ! mass tendency del4 + del4 frictional heating (dAH-dBH)
      70             :   real(r8) :: dMdt_residual             ! mass tendency residual (time truncation errors) 
      71             :   real(r8) :: dMdt_phys_total_in_dyn    ! mass tendency physics total in dycore
      72             :   real(r8) :: dMdt_PDC                  ! mass tendency physics-dynamics coupling
      73             :   !
      74             :   ! energy budgets dynamics
      75             :   !
      76             :   real(r8) :: dEdt_floating_dyn         ! dE/dt floating dynamics (dAL-dBL)
      77             :   real(r8) :: dEdt_vert_remap           ! dE/dt vertical remapping (dAR-dAD)
      78             :   real(r8) :: dEdt_del4                 ! dE/dt del4 (dCH-dBH)
      79             :   real(r8) :: dEdt_del4_fric_heat       ! dE/dt del4 frictional heating (dAH-dCH)
      80             :   real(r8) :: dEdt_del4_tot             ! dE/dt del4 + del4 fricitional heating (dAH-dBH)
      81             :   real(r8) :: dEdt_del2_sponge          ! dE/dt del2 sponge (dAS-dBS)
      82             :   real(r8) :: dEdt_del2_del4_tot        ! dE/dt explicit diffusion total
      83             :   real(r8) :: dEdt_residual             ! dE/dt residual (dEdt_floating_dyn-dEdt_del2_del4_tot)
      84             :   real(r8) :: dEdt_dycore_dyn           ! dE/dt adiabatic dynamical core (calculated in dycore)
      85             :   !
      86             :   ! physics-dynamics coupling variables
      87             :   !
      88             :   real(r8) :: E_dBF(4)                  ! E of dynamics state at the end of dycore integration (on dycore deomposition)
      89             :   real(r8) :: E_dyBF(4)                 ! E of physics state using dycore E
      90             : 
      91             : 
      92             :   real(r8) :: diff, tmp                 ! dummy variables
      93             :   integer  :: m_cnst, i
      94             :   character(LEN=*), parameter :: fmt  = "(a40,a15,a1,F6.2,a1,F6.2,a1,E10.2,a5)"
      95             :   character(LEN=*), parameter :: fmtf = "(a48,F8.4,a6)"
      96             :   character(LEN=*), parameter :: fmtm = "(a48,E8.2,a9)"
      97             :   character(LEN=15)           :: str(4)
      98             :   character(LEN=5)            :: pf     ! pass or fail identifier
      99             :   !--------------------------------------------------------------------------------------
     100             :   
     101      369408 :   if (masterproc .and. thermo_budget_history .and. hstwr(thermo_budget_histfile_num)) then
     102           0 :     idx(1) = teidx !total energy index
     103           0 :     idx(2) = seidx !enthaly index
     104           0 :     idx(3) = keidx !kinetic energy index
     105           0 :     idx(4) = poidx !surface potential energy index
     106           0 :     str(1) = "(total        )"
     107           0 :     str(2) = "(enthalpy     )"
     108           0 :     str(3) = "(kinetic      )"
     109           0 :     str(4) = "(srf potential)"
     110           0 :     do i=1,4
     111             :       !
     112             :       ! CAM physics energy tendencies
     113             :       !
     114           0 :       call cam_budget_get_global('phAP-phBP',idx(i),dEdt_param_physE(i))
     115             :       call cam_budget_get_global('phBP-phBF',idx(i),dEdt_efix_physE(i))
     116             :       call cam_budget_get_global('phAM-phAP',idx(i),dEdt_dme_adjust_physE(i))
     117             :       call cam_budget_get_global('phAP-phBF',idx(i),dEdt_param_efix_physE(i))
     118             :       !
     119             :       ! CAM physics energy tendencies using dycore energy formula scaling
     120             :       ! temperature tendencies for consistency with CAM physics
     121             :       !
     122             :       call cam_budget_get_global('dyAP-dyBP',idx(i),dEdt_param_dynE(i))
     123             :       call cam_budget_get_global('dyBP-dyBF',idx(i),dEdt_efix_dynE(i))
     124             :       call cam_budget_get_global('dyAM-dyAP',idx(i),dEdt_dme_adjust_dynE(i))
     125             :       call cam_budget_get_global('dyAP-dyBF',idx(i),dEdt_param_efix_dynE(i))
     126             :       call cam_budget_get_global('dyAM-dyBF',idx(i),dEdt_phys_total_dynE(i))
     127             :       call cam_budget_get_global('dyBF'     ,idx(i),E_dyBF(i))!state beginning physics
     128             :       !
     129             :       ! CAM physics energy tendencies in dynamical core
     130             :       !
     131             :       call cam_budget_get_global('dBD-dAF',idx(i),dEdt_phys_total_in_dyn(i))
     132           0 :       call cam_budget_get_global('dBF'    ,idx(i),E_dBF(i))  !state passed to physics
     133             :     end do
     134             : 
     135           0 :     call cam_budget_get_global('dAL-dBL',teidx,dEdt_floating_dyn)
     136           0 :     call cam_budget_get_global('dAR-dAD',teidx,dEdt_vert_remap)
     137           0 :     dEdt_dycore_dyn = dEdt_floating_dyn+dEdt_vert_remap
     138             : 
     139           0 :     call cam_budget_get_global('dCH-dBH',teidx,dEdt_del4)
     140           0 :     call cam_budget_get_global('dAH-dCH',teidx,dEdt_del4_fric_heat)
     141           0 :     call cam_budget_get_global('dAH-dBH',teidx,dEdt_del4_tot)
     142           0 :     call cam_budget_get_global('dAS-dBS',teidx,dEdt_del2_sponge)
     143           0 :     dEdt_del2_del4_tot      = dEdt_del4_tot+dEdt_del2_sponge
     144           0 :     dEdt_residual           = dEdt_floating_dyn-dEdt_del2_del4_tot
     145             : 
     146           0 :     write(iulog,*)" "
     147           0 :     write(iulog,*)"======================================================================"
     148           0 :     write(iulog,*)"Total energy diagnostics introduced in Lauritzen and Williamson (2019)"
     149           0 :     write(iulog,*)"(DOI:10.1029/2018MS001549)"
     150           0 :     write(iulog,*)"======================================================================"
     151           0 :     write(iulog,*)" "
     152           0 :     write(iulog,*)"Globally and vertically integrated total energy (E) diagnostics are"
     153           0 :     write(iulog,*)"computed at various points in the physics and dynamics loops to compute"
     154           0 :     write(iulog,*)"energy tendencies (dE/dt) and check for consistency (e.g., is E of"
     155           0 :     write(iulog,*)"state passed to physics computed using dycore state variables the same"
     156           0 :     write(iulog,*)"E of the state in the beginning of physics computed using the physics"
     157           0 :     write(iulog,*)"representation of the state)"
     158           0 :     write(iulog,*)" "
     159           0 :     write(iulog,*)"Energy stages in physics:"
     160           0 :     write(iulog,*)"-------------------------"
     161           0 :     write(iulog,*)" "
     162           0 :     write(iulog,*)"  xxBF: state passed to parameterizations, before energy fixer"
     163           0 :     write(iulog,*)"  xxBP: after energy fixer, before parameterizations"
     164           0 :     write(iulog,*)"  xxAP: after last phys_update in parameterizations and state "
     165           0 :     write(iulog,*)"        saved for energy fixer"
     166           0 :     write(iulog,*)"  xxAM: after dry mass adjustment"
     167           0 :     write(iulog,*)"  history files saved off here"
     168           0 :     write(iulog,*)" "
     169           0 :     write(iulog,*)"where xx='ph','dy' "
     170           0 :     write(iulog,*)" "
     171           0 :     write(iulog,*)"Suffix ph is CAM physics total energy"
     172           0 :     write(iulog,*)"(eq. 111 in Lauritzen et al. 2022; 10.1029/2022MS003117)"
     173           0 :     write(iulog,*)" "
     174           0 :     write(iulog,*)"Suffix dy is dycore energy computed in CAM physics using"
     175           0 :     write(iulog,*)"CAM physics state variables"
     176           0 :     write(iulog,*)" "
     177           0 :     write(iulog,*)" "
     178           0 :     write(iulog,*)"Energy stages in dynamics (specific to the SE dycore)"
     179           0 :     write(iulog,*)"-----------------------------------------------------"
     180           0 :     write(iulog,*)" "
     181           0 :     write(iulog,*)"suffix (d)"
     182           0 :     write(iulog,*)"dED: state from end of previous dynamics (= pBF + time sampling)"
     183           0 :     write(iulog,*)"   loop over vertical remapping and physics dribbling -------- (nsplit) -------"
     184           0 :     write(iulog,*)"            (dribbling and remapping always done together)                    |"
     185           0 :     write(iulog,*)"          dAF: state from previous remapping                                  |"
     186           0 :     write(iulog,*)"          dBD: state after physics dribble, before dynamics                   |"
     187           0 :     write(iulog,*)"          loop over vertical Lagrangian dynamics --------rsplit-------------  |"
     188           0 :     write(iulog,*)"              dynamics here                                                |  |"
     189           0 :     write(iulog,*)"              loop over hyperviscosity ----------hypervis_sub------------  |  |"
     190           0 :     write(iulog,*)"                 dBH   state before hyperviscosity                      |  |  |"
     191           0 :     write(iulog,*)"                 dCH   state after hyperviscosity                       |  |  |"
     192           0 :     write(iulog,*)"                 dAH   state after hyperviscosity momentum heating      |  |  |"
     193           0 :     write(iulog,*)"              end hyperviscosity loop -----------------------------------  |  |"
     194           0 :     write(iulog,*)"              dBS   state before del2 sponge                            |  |  |"
     195           0 :     write(iulog,*)"              dAS   state after del2+mom heating sponge                 |  |  |"
     196           0 :     write(iulog,*)"          end of vertical Lagrangian dynamics loop -------------------------  |"
     197           0 :     write(iulog,*)"      dAD  state after dynamics, before vertical remapping                    |"
     198           0 :     write(iulog,*)"      dAR     state after vertical remapping                                  |"
     199           0 :     write(iulog,*)"   end of remapping loop ------------------------------------------------------"
     200           0 :     write(iulog,*)"dBF  state passed to parameterizations = state after last remapping            "
     201           0 :     write(iulog,*)" "
     202           0 :     write(iulog,*)" "
     203           0 :     write(iulog,*)"FYI: all difference (diff) below are absolute normalized differences"
     204           0 :     write(iulog,*)" "
     205           0 :     write(iulog,*)"Consistency check 0:"
     206           0 :     write(iulog,*)"--------------------"
     207           0 :     write(iulog,*)" "
     208           0 :     write(iulog,*)"For energetic consistency we require that dE/dt [W/m^2] from energy "
     209           0 :     write(iulog,*)"fixer and all parameterizations computed using physics E and"
     210           0 :     write(iulog,*)"dycore in physics E are the same! Checking:"
     211           0 :     write(iulog,*)" "
     212           0 :     write(iulog,*)  "                                                        xx=ph   xx=dy  norm. diff."
     213           0 :     write(iulog,*)  "                                                        -----   -----  -----------"
     214           0 :     do i=1,4
     215           0 :       diff = abs_diff(dEdt_efix_physE(i),dEdt_efix_dynE(i),pf=pf)
     216           0 :       write(iulog,fmt)"dE/dt energy fixer          (xxBP-xxBF) ",str(i)," ",dEdt_efix_physE(i), " ", &
     217           0 :                        dEdt_efix_dynE(i)," ",diff,pf
     218           0 :       diff = abs_diff(dEdt_param_physE(i),dEdt_param_dynE(i),pf=pf)
     219           0 :       write(iulog,fmt)"dE/dt all parameterizations (xxAP-xxBP) ",str(i)," ",dEdt_param_physE(i)," ", &
     220           0 :                       dEdt_param_dynE(i)," ",diff,pf
     221           0 :       write(iulog,*) " "
     222           0 :       if (diff>eps) then
     223           0 :         write(iulog,*)"FAIL"
     224           0 :         call endrun(subname//"dE/dt's in physics inconsistent")
     225             :       end if
     226             :     end do
     227           0 :     write(iulog,*)" "
     228           0 :     write(iulog,*)" "
     229           0 :     write(iulog,*)"dE/dt from dry-mass adjustment will differ if dynamics and physics use"
     230           0 :     write(iulog,*)"different energy definitions! Checking:"
     231           0 :     write(iulog,*)" "
     232           0 :     write(iulog,*)  "                                                        xx=ph   xx=dy  diff"
     233           0 :     write(iulog,*)  "                                                        -----   -----  ----"
     234           0 :     do i=1,4
     235           0 :       diff = dEdt_dme_adjust_physE(i)-dEdt_dme_adjust_dynE(i)
     236           0 :       write(iulog,fmt)"dE/dt dry mass adjustment   (xxAM-xxAP) ",str(i)," ",dEdt_dme_adjust_physE(i)," ", &
     237           0 :                       dEdt_dme_adjust_dynE(i)," ",diff
     238             :     end do
     239           0 :     write(iulog,*)" "
     240           0 :     write(iulog,*)" "
     241             :     !
     242             :     ! these diagnostics only make sense time-step to time-step
     243             :     !
     244           0 :     write(iulog,*)" "
     245           0 :     write(iulog,*)"Some energy budget observations:"
     246           0 :     write(iulog,*)"--------------------------------"
     247           0 :     write(iulog,*)" "
     248           0 :     write(iulog,*)"Note that total energy fixer fixes:"
     249           0 :     write(iulog,*) " "
     250           0 :     write(iulog,*) "-dE/dt energy fixer(t=n) = dE/dt dry mass adjustment (t=n-1) +"
     251           0 :     write(iulog,*) "                      dE/dt adiabatic dycore (t=n-1)         +"
     252           0 :     write(iulog,*) "                      dE/dt physics-dynamics coupling errors (t=n-1)"
     253           0 :     write(iulog,*) " "
     254           0 :     write(iulog,*) "(equation 23 in Lauritzen and Williamson (2019))"
     255           0 :     write(iulog,*) " "
     256             : 
     257           0 :     tmp = previous_dEdt_phys_dyn_coupl_err+previous_dEdt_adiabatic_dycore+previous_dEdt_dry_mass_adjust
     258           0 :     diff = abs_diff(-dEdt_efix_dynE(1),tmp,pf)
     259           0 :     if (.not.use_cslam) then
     260           0 :       write(iulog,*) "Check if that is the case:", pf, diff
     261           0 :       write(iulog,*) " "
     262           0 :       if (abs(diff)>eps) then
     263           0 :         write(iulog,*) "dE/dt energy fixer(t=n)                        = ",dEdt_efix_dynE(1)
     264           0 :         write(iulog,*) "dE/dt dry mass adjustment (t=n-1)              = ",previous_dEdt_dry_mass_adjust
     265           0 :         write(iulog,*) "dE/dt adiabatic dycore (t=n-1)                 = ",previous_dEdt_adiabatic_dycore
     266           0 :         write(iulog,*) "dE/dt physics-dynamics coupling errors (t=n-1) = ",previous_dEdt_phys_dyn_coupl_err
     267             :       end if
     268             :     else
     269           0 :       previous_dEdt_phys_dyn_coupl_err = dEdt_efix_dynE(1)+previous_dEdt_dry_mass_adjust+previous_dEdt_adiabatic_dycore
     270           0 :       write(iulog,*) "dE/dt energy fixer(t=n)                        = ",dEdt_efix_dynE(1)
     271           0 :       write(iulog,*) "dE/dt dry mass adjustment (t=n-1)              = ",previous_dEdt_dry_mass_adjust
     272           0 :       write(iulog,*) "dE/dt adiabatic dycore (t=n-1)                 = ",previous_dEdt_adiabatic_dycore
     273           0 :       write(iulog,*) "dE/dt physics-dynamics coupling errors (t=n-1) = ",previous_dEdt_phys_dyn_coupl_err
     274           0 :       write(iulog,*) " "
     275           0 :       write(iulog,*) "Note: when running CSLAM the physics-dynamics coupling error is diagnosed"
     276           0 :       write(iulog,*) "      (using equation above) rather than explicitly computed"
     277           0 :       write(iulog,*) " "
     278           0 :       write(iulog,*) " "
     279           0 :       write(iulog,*) "Physics-dynamics coupling errors include: "
     280           0 :       write(iulog,*) " "
     281           0 :       write(iulog,*) " -dE/dt adiabatic dycore is computed on GLL grid;"
     282           0 :       write(iulog,*) " error in mapping to physics grid"
     283           0 :       write(iulog,*) " -dE/dt physics tendencies mapped to GLL grid"
     284           0 :       write(iulog,*) " (tracer tendencies mapped non-conservatively!)"
     285           0 :       write(iulog,*) " -dE/dt dynamics state mapped to GLL grid"
     286             :     end if
     287           0 :     write(iulog,*) ""
     288           0 :     if (.not.use_cslam) then
     289           0 :       dEdt_dycore_phys = -dEdt_efix_dynE(1)-previous_dEdt_phys_dyn_coupl_err-previous_dEdt_dry_mass_adjust
     290           0 :       write(iulog,*)               "Hence the dycore E dissipation estimated from energy fixer "
     291           0 :       write(iulog,'(A39,F6.2,A6)') "based on previous time-step values is ",dEdt_dycore_phys," W/M^2"
     292           0 :       write(iulog,*) " "
     293             :     end if
     294           0 :     write(iulog,*) " "
     295           0 :     write(iulog,*) "-------------------------------------------------------------------"
     296           0 :     write(iulog,*) " Consistency check 1: state passed to physics same as end dynamics?"
     297           0 :     write(iulog,*) "-------------------------------------------------------------------"
     298           0 :     write(iulog,*) " "
     299           0 :     write(iulog,*) "Is globally integrated total energy of state at the end of dynamics (dBF)"
     300           0 :     write(iulog,*) "and beginning of physics (using dynamics in physics energy; dyBF) the same?"
     301           0 :     write(iulog,*) ""
     302           0 :     if (.not.use_cslam) then
     303           0 :       if (abs(E_dyBF(1))>eps) then
     304           0 :         diff = abs_diff(E_dBF(1),E_dyBF(1))
     305           0 :         if (abs(diff)<eps) then
     306           0 :           write(iulog,'(A23,E8.3)')"yes. (dBF-dyBF)/dyBF =",diff
     307             :         else
     308           0 :           write(iulog,*)"Error in physics dynamics coupling!"
     309           0 :           write(iulog,*)" "
     310           0 :           do i=1,4
     311           0 :             write(iulog,*) str(i),":"
     312           0 :             write(iulog,*) "======"
     313           0 :             diff = abs_diff(E_dBF(i),E_dyBF(i),pf=pf)
     314           0 :             write(iulog,*) "diff, E_dBF, E_dyBF ",diff,E_dBF(i),E_dyBF(i)
     315           0 :             write(iulog,*) " "
     316             :           end do
     317           0 :           call endrun(subname//"Error in physics dynamics coupling")
     318             :         end if
     319             :       end if
     320             :     else
     321           0 :       write(iulog,*)" "
     322           0 :       write(iulog,*)"Since you are using a separate physics grid, the state in dynamics"
     323           0 :       write(iulog,*)"will not be the same on the physics grid since it is"
     324           0 :       write(iulog,*)"interpolated from the dynamics to the physics grid"
     325           0 :       write(iulog,*)" "
     326           0 :       do i=1,4
     327           0 :         write(iulog,*) str(i),":"
     328           0 :         write(iulog,*) "======"
     329           0 :         diff = abs_diff(E_dBF(i),E_dyBF(i),pf=pf)
     330           0 :         write(iulog,*) "diff, E_dBF, E_dyBF ",diff,E_dBF(i),E_dyBF(i)
     331           0 :         write(iulog,*) " "
     332             :       end do
     333             :     end if
     334             : 
     335           0 :     write(iulog,*)" "
     336           0 :     write(iulog,*)"-------------------------------------------------------------------------"
     337           0 :     write(iulog,*)" Consistency check 2: total energy increment in dynamics same as physics?"
     338           0 :     write(iulog,*)"-------------------------------------------------------------------------"
     339           0 :     write(iulog,*)" "
     340           0 :     if (.not.use_cslam) then
     341           0 :       previous_dEdt_phys_dyn_coupl_err = dEdt_phys_total_in_dyn(1)-dEdt_phys_total_dynE(1)
     342           0 :       diff = abs_diff(dEdt_phys_total_dynE(1),dEdt_phys_total_in_dyn(1),pf=pf)
     343           0 :       write(iulog,'(A40,E8.2,A7,A5)')" dE/dt physics-dynamics coupling errors       ",diff," W/M^2 ",pf
     344           0 :       if (abs(diff)>eps) then
     345             :         !
     346             :         ! if errors print details
     347             :         !
     348           0 :         if (ftype==1) then
     349           0 :           write(iulog,*) ""
     350           0 :           write(iulog,*) " You are using ftype==1 so physics-dynamics coupling errors should be round-off!"
     351           0 :           write(iulog,*) ""
     352           0 :           write(iulog,*) " Because of failure provide detailed diagnostics below:"
     353           0 :           write(iulog,*) ""
     354             :         else
     355           0 :           write(iulog,*) ""
     356           0 :           write(iulog,*) " Since ftype<>1 there are physics dynamics coupling errors"
     357           0 :           write(iulog,*) ""
     358           0 :           write(iulog,*) " Break-down below:"
     359           0 :           write(iulog,*) ""
     360             :         end if
     361             :         
     362           0 :         do i=1,4
     363           0 :           write(iulog,*) str(i),":"
     364           0 :           write(iulog,*) "======"
     365           0 :           diff = abs_diff(dEdt_phys_total_dynE(i),dEdt_phys_total_in_dyn(i),pf=pf)
     366           0 :           write(iulog,*) "dE/dt physics-dynamics coupling errors (diff) ",diff
     367           0 :           write(iulog,*) "dE/dt physics total in dynamics (dBD-dAF)     ",dEdt_phys_total_in_dyn(i)
     368           0 :           write(iulog,*) "dE/dt physics total in physics  (dyAM-dyBF)   ",dEdt_phys_total_dynE(i)
     369           0 :           write(iulog,*) " "
     370           0 :           write(iulog,*) "      physics total = parameterizations + efix + dry-mass adjustment"
     371           0 :           write(iulog,*) " "
     372             :         end do
     373             : !   Temporarily disable endrun until energy bias for consistancy check 2 is better understood.
     374             : !        if (ftype==1) then
     375             : !          call endrun(subname//"Physics-dynamics coupling error. See atm.log")
     376             : !        end if
     377             :       end if
     378             :     else
     379           0 :       write(iulog,'(a47,F6.2,a6)')" dE/dt physics tendency in dynamics (dBD-dAF)   ",dEdt_phys_total_in_dyn(1)," W/M^2"
     380           0 :       write(iulog,'(a47,F6.2,a6)')" dE/dt physics tendency in physics  (dyAM-dyBF) ",dEdt_phys_total_dynE(1)," W/M^2"
     381           0 :       write(iulog,*)" "
     382           0 :       write(iulog,*) " When runnig with a physics grid this consistency check does not make sense"
     383           0 :       write(iulog,*) " since it is computed on the GLL grid whereas we enforce energy conservation"
     384           0 :       write(iulog,*) " on the physics grid. To assess the errors of running dynamics on GLL"
     385           0 :       write(iulog,*) " grid, tracers on CSLAM grid and physics on physics grid we use the energy"
     386           0 :       write(iulog,*) " fixer check from above:"
     387           0 :       write(iulog,*) " "
     388           0 :       write(iulog,*) " dE/dt physics-dynamics coupling errors (t=n-1) =",previous_dEdt_phys_dyn_coupl_err
     389           0 :       write(iulog,*) ""
     390             :     end if
     391           0 :     write(iulog,*)" "
     392           0 :     write(iulog,*)"------------------------------------------------------------"
     393           0 :     write(iulog,*)" SE dycore energy tendencies"
     394           0 :     write(iulog,*)"------------------------------------------------------------"
     395           0 :     write(iulog,*)" "
     396           0 :     write(iulog,fmtf)"   dE/dt dycore                                  ",dEdt_dycore_dyn," W/M^2"
     397           0 :     write(iulog,*)" "
     398           0 :     write(iulog,*)"Adiabatic dynamics can be divided into quasi-horizontal and vertical remapping: "
     399           0 :     write(iulog,*)" "
     400           0 :     write(iulog,fmtf)"   dE/dt floating dynamics           (dAD-dBD)   ",dEdt_floating_dyn," W/M^2"
     401           0 :     write(iulog,fmtf)"   dE/dt vertical remapping          (dAR-dAD)   ",dEdt_vert_remap," W/M^2"
     402             : 
     403           0 :     write(iulog,*) " "
     404           0 :     write(iulog,*) "Breakdown of floating dynamics:"
     405           0 :     write(iulog,*) " "
     406           0 :     write(iulog,fmtf)"   dE/dt hypervis del4               (dCH-dBH)   ",dEdt_del4,          " W/M^2"
     407           0 :     write(iulog,fmtf)"   dE/dt hypervis frictional heating (dAH-dCH)   ",dEdt_del4_fric_heat," W/M^2"
     408           0 :     write(iulog,fmtf)"   dE/dt hypervis del4 total         (dAH-dBH)   ",dEdt_del4_tot, " W/M^2"
     409           0 :     write(iulog,fmtf)"   dE/dt hypervis sponge del2        (dAS-dBS)   ",dEdt_del2_sponge,   " W/M^2"
     410           0 :     write(iulog,fmtf)"   dE/dt explicit diffusion total                ",dEdt_del2_del4_tot,    " W/M^2"
     411           0 :     write(iulog,*) " "
     412           0 :     write(iulog,fmtf)"   dE/dt residual (time-truncation errors,...)   ",dEdt_residual,      " W/M^2"
     413           0 :     write(iulog,*)" "
     414           0 :     write(iulog,*)" "
     415           0 :     write(iulog,*)"------------------------------------------------------------"
     416           0 :     write(iulog,*)"Tracer mass budgets"
     417           0 :     write(iulog,*)"------------------------------------------------------------"
     418           0 :     write(iulog,*)" "
     419           0 :     write(iulog,*)"Below the physics-dynamics coupling error is computed as    "
     420           0 :     write(iulog,*)"dMASS/dt physics tendency in dycore (dBD-dAF) minus"
     421           0 :     write(iulog,*)"dMASS/dt total physics              (pAM-pBF)"
     422           0 :     write(iulog,*)" "
     423           0 :     write(iulog,*)" "
     424           0 :     do m_cnst=1,thermo_budget_num_vars
     425           0 :       if (thermo_budget_vars_massv(m_cnst)) then
     426           0 :         write(iulog,*)thermo_budget_vars_descriptor(m_cnst)
     427           0 :         write(iulog,*)"------------------------------"
     428           0 :         call cam_budget_get_global('phBP-phBF',m_cnst,dMdt_efix)
     429           0 :         call cam_budget_get_global('phAM-phAP',m_cnst,dMdt_dme_adjust)
     430           0 :         call cam_budget_get_global('phAP-phBP',m_cnst,dMdt_parameterizations)
     431           0 :         call cam_budget_get_global('phAM-phBF',m_cnst,dMdt_phys_total)
     432             :         !
     433             :         ! total energy fixer should not affect mass - checking
     434             :         !
     435           0 :         if (abs(dMdt_efix)>eps_mass) then
     436           0 :           write(iulog,*) "dMASS/dt energy fixer        (pBP-pBF)           ",dMdt_efix," Pa/m^2/s"
     437           0 :           write(iulog,*) "ERROR: Mass not conserved in energy fixer. ABORT"      
     438           0 :           call endrun(subname//"Mass not conserved in energy fixer. See atm.log")
     439             :         endif
     440             :         !
     441             :         ! dry-mass adjustmnt should not affect mass - checking
     442             :         !
     443           0 :         if (abs(dMdt_dme_adjust)>eps_mass) then
     444           0 :           write(iulog,*)"dMASS/dt dry mass adjustment (pAM-pAP) ",dMdt_dme_adjust," Pa/m^2/s"
     445           0 :           write(iulog,*) "ERROR: Mass not conserved in dry mass adjustment. ABORT"
     446           0 :           call endrun(subname//"Mass not conserved in dry mass adjustment. See atm.log")
     447             :         end if
     448             :         !
     449             :         ! all of the mass-tendency should come from parameterization - checking
     450             :         !
     451           0 :         if (abs(dMdt_parameterizations-dMdt_phys_total)>eps_mass) then
     452           0 :           write(iulog,*) "Error: dMASS/dt parameterizations (pAP-pBP) .ne. dMASS/dt physics total (pAM-pBF)"
     453           0 :           write(iulog,*) "dMASS/dt parameterizations   (pAP-pBP) ",dMdt_parameterizations," Pa/m^2/s"
     454           0 :           write(iulog,*) "dMASS/dt physics total       (pAM-pBF) ",dMdt_phys_total," Pa/m^2/s"
     455           0 :           call endrun(subname//"mass change not only due to parameterizations. See atm.log")
     456             :         end if
     457           0 :         write(iulog,*)"  "
     458             :         !
     459             :         ! detailed mass budget in dynamical core
     460             :         !
     461           0 :         if (is_cam_budget('dAD').and.is_cam_budget('dBD').and.is_cam_budget('dAR').and.is_cam_budget('dCH')) then
     462           0 :           call cam_budget_get_global('dAL-dBL',m_cnst,dMdt_floating_dyn)
     463           0 :           call cam_budget_get_global('dAR-dAD',m_cnst,dMdt_vert_remap)
     464           0 :           tmp  = dMdt_floating_dyn+dMdt_vert_remap
     465           0 :           diff = abs_diff(tmp,0.0_r8,pf=pf)
     466           0 :           write(iulog,fmtm)"   dMASS/dt total adiabatic dynamics             ",diff,pf
     467             :           !
     468             :           ! check for mass-conservation in the adiabatic dynamical core -
     469             :           ! if not conserved provide detailed break-down
     470             :           !
     471           0 :           if (abs(diff)>eps_mass) then
     472           0 :             write(iulog,*) "Error: mass non-conservation in dynamical core"
     473           0 :             write(iulog,*) "(detailed budget below)"
     474           0 :             write(iulog,*) " "
     475           0 :             write(iulog,*)"dMASS/dt 2D dynamics            (dAL-dBL) ",dMdt_floating_dyn," Pa/m^2/s"
     476           0 :             write(iulog,*)"dE/dt vertical remapping        (dAR-dAD) ",dMdt_vert_remap
     477           0 :             write(iulog,*)" "
     478           0 :             write(iulog,*)"Breakdown of 2D dynamics:"
     479           0 :             write(iulog,*)" "
     480           0 :             call cam_budget_get_global('dAH-dCH',m_cnst,dMdt_del4_fric_heat)
     481           0 :             call cam_budget_get_global('dAH-dBH',m_cnst,dMdt_del4_tot)
     482           0 :             write(iulog,*)"dMASS/dt hypervis               (dAH-dBH) ",dMdt_del4_tot," Pa/m^2/s"
     483           0 :             write(iulog,*)"dMASS/dt frictional heating     (dAH-dCH) ",dMdt_del4_fric_heat," Pa/m^2/s"
     484           0 :             dMdt_residual = dMdt_floating_dyn-dMdt_del4_tot
     485           0 :             write(iulog,*)"dMASS/dt residual (time truncation errors)",dMdt_residual," Pa/m^2/s"
     486             :           end if
     487             :         end if
     488           0 :         if (is_cam_budget('dBD').and.is_cam_budget('dAF')) then
     489             :           !
     490             :           ! check if mass change in physics is the same as dynamical core
     491             :           !
     492           0 :           call cam_budget_get_global('dBD-dAF',m_cnst,dMdt_phys_total_in_dyn)
     493           0 :           dMdt_PDC = dMdt_phys_total-dMdt_phys_total_in_dyn
     494           0 :           write(iulog,fmtm)"   Mass physics-dynamics coupling error          ",dMdt_PDC," Pa/m^2/s"
     495           0 :           write(iulog,*)" "
     496           0 :           if (abs(dMdt_PDC)>eps_mass) then
     497           0 :             write(iulog,fmtm)"   dMASS/dt physics tendency in dycore (dBD-dAF) ",dMdt_phys_total_in_dyn," Pa/m^2/s"
     498           0 :             write(iulog,fmtm)"   dMASS/dt total physics                        ",dMdt_phys_total," Pa/m^2/s"
     499             :           end if
     500             :         end if
     501             :       end if
     502             :     end do
     503             :     !
     504             :     ! save adiabatic dycore dE/dt and dry-mass adjustment to avoid samping error
     505             :     !
     506           0 :     previous_dEdt_adiabatic_dycore = dEdt_dycore_dyn
     507           0 :     previous_dEdt_dry_mass_adjust  = dEdt_dme_adjust_dynE(1)
     508             :   end if
     509      369408 : end subroutine print_budget
     510             : !=========================================================================================
     511           0 : function abs_diff(a,b,pf)
     512             :   real(r8),                   intent(in) :: a,b
     513             :   character(LEN=5), optional, intent(out):: pf
     514             :   real(r8)                               :: abs_diff
     515           0 :   if (abs(b)>eps) then
     516           0 :     abs_diff = abs((b-a)/b)
     517             :   else
     518           0 :     abs_diff = abs(b-a)
     519             :   end if
     520           0 :   If (present(pf)) then
     521           0 :     if (abs_diff>eps) then
     522           0 :       pf = ' FAIL'
     523             :     else
     524           0 :       pf = ' PASS'
     525             :     end if
     526             :   end if
     527      369408 : end function abs_diff
     528             : end module dycore_budget

Generated by: LCOV version 1.14