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 22:39:59 |  | 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