From 09720b32684ef24cccd389fe39cf3bd6f9f32527 Mon Sep 17 00:00:00 2001 From: sit23 Date: Tue, 11 Aug 2026 11:55:56 +0100 Subject: [PATCH 1/5] Merge orographic (mg_drag) and non-orographic (cg_drag) GWD scheme changes from Ross Castle's fork Brings in Ross's PhD-era physics changes to the mountain gravity wave drag scheme (from ross/realistic_mgdrag) and the convective gravity wave drag scheme (from ross/gwdp_chaim_steven), taking mg_drag.f90 and cg_drag.f90 wholesale from their respective branches. Both branches independently added is/ie/js/je domain-decomposition bounds arguments to mg_drag_init/cg_drag_init (needed for array allocation in mg_drag_init). damping_driver.f90 and its damping_driver_init call site in idealized_moist_phys.F90 are hand-reconciled to thread these through to both schemes, since master has evolved independently from either branch's fork point and a naive merge would have conflicted or silently dropped one side's wiring. Not included: cg_drag_current.f90 (an unwired backup of the pre-edit file, per its own header) and an unrelated fftw.F90 path_names addition from the mg_drag branch. --- src/atmos_param/cg_drag/cg_drag.f90 | 642 ++++++++++-------- .../damping_driver/damping_driver.f90 | 31 +- src/atmos_param/mg_drag/mg_drag.f90 | 247 +++++-- .../driver/solo/idealized_moist_phys.F90 | 2 +- 4 files changed, 576 insertions(+), 346 deletions(-) diff --git a/src/atmos_param/cg_drag/cg_drag.f90 b/src/atmos_param/cg_drag/cg_drag.f90 index 41bee433c..879021181 100644 --- a/src/atmos_param/cg_drag/cg_drag.f90 +++ b/src/atmos_param/cg_drag/cg_drag.f90 @@ -1,3 +1,4 @@ +!Chaim cg_drag with Stepehens cg_drag edits module cg_drag_mod use fms_mod, only: fms_init, mpp_pe, mpp_root_pe, & @@ -29,11 +30,11 @@ module cg_drag_mod private !--------------------------------------------------------------------- -! cg_drag_mod computes the convective gravity wave forcing on -! the zonal flow. the parameterization is described in Alexander and -! Dunkerton [JAS, 15 December 1999]. +! cg_drag_mod computes the convective gravity wave forcing on +! the zonal flow. the parameterization is described in Alexander and +! Dunkerton [JAS, 15 December 1999]. !-------------------------------------------------------------------- - + !--------------------------------------------------------------------- !----------- ****** VERSION NUMBER ******* --------------------------- @@ -63,7 +64,7 @@ module cg_drag_mod !!$integer :: vers, old_time_step !wfc++ Addition for regular use - integer, allocatable, dimension(:,:) :: source_level + integer, allocatable, dimension(:,:) :: source_level, damp_level real, allocatable, dimension(:,:) :: source_amp real, allocatable, dimension(:,:,:) :: gwd_u, gwd_v @@ -79,68 +80,100 @@ module cg_drag_mod ! at 00Z and calculations are not done ! every time step -real :: source_level_pressure= 315.e+02 - ! highest model level with pressure +real :: source_level_pressure= 315.e+02 + ! highest model level with pressure ! greater than this value (or sigma ! greater than this value normalized ! by 1013.25 hPa) will be the gravity - ! wave source level at the equator + ! wave source level at the equator ! [ Pa ] -integer :: nk=1 ! number of wavelengths contained in +real :: damp_level_pressure=0.8e+02 + ! added by cig, feb 27, 2017. any waves reaching the top level will be deposited down to this level +integer :: nk=1 ! number of wavelengths contained in ! the gravity wave spectrum real :: cmax=99.6 ! maximum phase speed in gravity wave ! spectrum [ m/s ] -real :: dc=1.2 ! gravity wave spectral resolution +real :: dc=1.2 ! gravity wave spectral resolution ! [ m/s ] ! previous values: 0.6 -real :: Bt_0=.004 ! sum across the wave spectrum of +real :: Bt_0=.004 ! sum across the wave spectrum of ! the magnitude of total GW stress [Pa] - -real :: Bt_aug=.000 ! magnitude of momentum flux divided by density + +real :: Bt_aug=.000 ! magnitude of momentum flux divided by density real :: Bt_nh=.001 ! additional momentum stress for NH [Pa] real :: Bt_sh=-.001 ! additional momentum stress for SH [Pa] -real :: Bt_eq=.000 ! additional momentum stress at equator - CURRENTLY NOT USED! +! epg - 30.6.16 - I shifted these spectral parameters to the name list + +!--------------------------------------------------------------------- +!--------------------------------------------------------------------- +! wave spectrum parameters. +!--------------------------------------------------------------------- + +integer :: flag = 1 ! flag = 1 for peak flux at c = 0 + ! flag = 0 for peak flux at (c-u) = 0 +real :: Bw = 0.4 ! amplitude for the wide spectrum [ m^2/s^2 ] + ! ~ u'w' +real :: Bn = 0.0 ! amplitude for the narrow spectrum [ m^2/s^2 ] + ! ~ u'w'; previous values: 5.4 +real :: cw = 40.0 ! half-width for the wide c spectrum [ m/s ] + ! previous values: 50.0, 25.0 +real :: cwtropics = 40.0 ! half-width for the wide c spectrum [ m/s ] + ! previous values: 50.0, 25.0 +real :: cn = 2.0 ! half-width for the narrow c spectrum [ m/s ] + +real :: Bt_eq=.000 ! additional momentum stress at equator - CURRENTLY NOT USED! real :: Bt_eq_width=4.0 ! scaling for width of equtorial momentum flux (equator) CURRENTLY NOT USED! real :: phi0n = 30., phi0s = -30., dphin = 5., dphis = -5. -logical :: calculate_ked=.false. +!add by chaim jan 2017 +real :: weightminus2=0. + +real :: weightminus1=0. + +real :: weighttop=1. + +real :: kelvin_kludge=1. + +logical :: calculate_ked=.false. ! calculate ked diagnostic ? integer :: num_diag_pts_ij=0 ! number of diagnostic columns specif- ! ied by global (i,j) coordinates -integer :: num_diag_pts_latlon=0 +integer :: num_diag_pts_latlon=0 ! number of diagnostic columns ! specified by lat-lon coordinates integer, parameter :: MAX_PTS= 20 ! maximum number of diagnostic columns -integer, dimension(MAX_PTS) :: i_coords_gl=-100 - ! global i coordinates for ij - ! diagnostic columns -integer, dimension(MAX_PTS) :: j_coords_gl=-100 - ! global j coordinates for ij - ! diagnostic columns -real, dimension(MAX_PTS) :: lat_coords_gl=-999. - ! latitudes for latlon diagnostic +integer, dimension(MAX_PTS) :: i_coords_gl=-100 + ! global i coordinates for ij + ! diagnostic columns +integer, dimension(MAX_PTS) :: j_coords_gl=-100 + ! global j coordinates for ij + ! diagnostic columns +real, dimension(MAX_PTS) :: lat_coords_gl=-999. + ! latitudes for latlon diagnostic ! columns [degrees, -90. -> 90. ] -real, dimension(MAX_PTS) :: lon_coords_gl=-999. - ! longitudes for latlon diagnostic +real, dimension(MAX_PTS) :: lon_coords_gl=-999. + ! longitudes for latlon diagnostic ! columns [ degrees, 0. -> 360. ] namelist / cg_drag_nml / & cg_drag_freq, cg_drag_offset, & - source_level_pressure, & + source_level_pressure, damp_level_pressure, & nk, cmax, dc, Bt_0, Bt_aug, & Bt_sh, Bt_nh, Bt_eq, Bt_eq_width, & calculate_ked, & num_diag_pts_ij, num_diag_pts_latlon, & i_coords_gl, j_coords_gl, & lat_coords_gl, lon_coords_gl, & - phi0n,phi0s,dphin,dphis + phi0n,phi0s,dphin,dphis, Bw, Bn, cw, cwtropics, cn, flag, & + weightminus2, weightminus1, weighttop,kelvin_kludge + !-------------------------------------------------------------------- !-------- public data ----- @@ -154,16 +187,16 @@ module cg_drag_mod !-------------------------------------------------------------------- !mj remove restart stuff integer, dimension(3) :: restart_versions = (/ 1, 2, 3 /) ! v1 : -! v2 : +! v2 : ! v3 : Now use NetCDF for restart file. ! !-------------------------------------------------------------------- ! these arrays must be preserved across timesteps in case the ! parameterization is not called every timestep: ! -! gwd time tendency for u eqn due to gravity wave forcing +! gwd time tendency for u eqn due to gravity wave forcing ! [ m/s^2 ] -! ked effective eddy diffusion coefficient resulting from +! ked effective eddy diffusion coefficient resulting from ! gravity wave forcing [ m^2/s ] ! !-------------------------------------------------------------------- @@ -180,28 +213,18 @@ module cg_drag_mod !------------------------------------------------------------------- real, dimension(:), allocatable :: c0, kwv, k2 - -!--------------------------------------------------------------------- -! wave spectrum parameters. -!--------------------------------------------------------------------- integer :: nc ! number of wave speeds in spectrum ! (symmetric around c = 0) -integer :: flag = 1 ! flag = 1 for peak flux at c = 0 - ! flag = 0 for peak flux at (c-u) = 0 -real :: Bw = 0.4 ! amplitude for the wide spectrum [ m^2/s^2 ] - ! ~ u'w' -real :: Bn = 0.0 ! amplitude for the narrow spectrum [ m^2/s^2 ] - ! ~ u'w'; previous values: 5.4 -real :: cw = 40.0 ! half-width for the wide c spectrum [ m/s ] - ! previous values: 50.0, 25.0 -real :: cn = 2.0 ! half-width for the narrow c spectrum [ m/s ] -integer :: klevel_of_source +integer :: klevel_of_source, klevel_of_damp ! k index of the gravity wave source level at ! the equator in a standard atmosphere + ! also k index of level up to where mesosphere drag is dumped (cig, feb 27 2017) + + !--------------------------------------------------------------------- ! variables which control module calculations: -! +! ! cgdrag_alarm time remaining until next cg_drag calculation [ s ] ! !--------------------------------------------------------------------- @@ -211,29 +234,30 @@ module cg_drag_mod ! variables used with column diagnostics: ! ! diag_units output unit numbers -! num_diag_pts number of columns where diagnostics are desired +! num_diag_pts number of columns where diagnostics are desired ! column_diagnostics_desired ! column diagnostics are desired ? -! do_column_diagnostics -! a diagnostic column is in this jrow ? +! do_column_diagnostics +! a diagnostic column is in this jrow ? ! diag_lon longitude of diagnostic columns [ degrees ] ! diag_lat latiude of diagnostic columns [ degrees ] ! diag_i processor-based i index of diagnostic columns ! diag_j processor-based j index of diagnostic columns ! !-------------------------------------------------------------------- -integer :: num_diag_pts = 0 +integer :: num_diag_pts = 0 logical :: column_diagnostics_desired=.false. -integer, dimension(:), allocatable :: diag_units +integer, dimension(:), allocatable :: diag_units logical, dimension(:), allocatable :: do_column_diagnostics real, dimension(:), allocatable :: diag_lon, diag_lat -integer, dimension(:), allocatable :: diag_j, diag_i +integer, dimension(:), allocatable :: diag_j, diag_i !--------------------------------------------------------------------- ! variables for netcdf diagnostic fields. !--------------------------------------------------------------------- integer :: id_kedx_cgwd, id_kedy_cgwd, id_bf_cgwd, & - id_gwfx_cgwd, id_gwfy_cgwd +! id_gwfx_cgwd, id_gwfy_cgwd + id_gwfx_cgwd, id_gwfy_cgwd, id_source_level real :: missing_value = -999. character(len=7) :: mod_name = 'cg_drag' @@ -256,8 +280,8 @@ module cg_drag_mod !#################################################################### -subroutine cg_drag_init (lonb, latb, pref, Time, axes) - +!subroutine cg_drag_init (lonb, latb, pref, Time, axes) RC 5/12/19 +subroutine cg_drag_init (is, ie, js, je, lonb, latb, pref, Time, axes) !------------------------------------------------------------------- ! cg_drag_init is the constructor for cg_drag_mod. @@ -268,6 +292,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) real, dimension(:), intent(in) :: lonb, latb, pref integer, dimension(4), intent(in) :: axes type(time_type), intent(in) :: Time +integer, intent(in) :: is, ie, js, je !RC 5/12/19 !------------------------------------------------------------------- !------------------------------------------------------------------- @@ -284,7 +309,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) !------------------------------------------------------------------ !------------------------------------------------------------------- -! local variables: +! local variables: integer :: unit, ierr, io, logunit integer :: n, i, j, k @@ -297,11 +322,11 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) ! real :: lat(size(lonb,1) - 1, size(latb,2) - 1) real :: lat(size(lonb) - 1, size(latb) - 1) !------------------------------------------------------------------- -! local variables: -! -! unit unit number for nml file -! ierr error return flag -! io error return code +! local variables: +! +! unit unit number for nml file +! ierr error return flag +! io error return code ! n loop index ! k loop index ! idf number of i points on this processor @@ -323,7 +348,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) call diag_manager_init call constants_init #ifdef COL_DIAG - call column_diagnostics_init +! call column_diagnostics_init #endif SKIP !--------------------------------------------------------------------- ! read namelist. @@ -345,7 +370,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) if (mpp_pe() == mpp_root_pe()) write (logunit, nml=cg_drag_nml) !------------------------------------------------------------------- -! define the grid dimensions. idf and jdf are the (i,j) dimensions of +! define the grid dimensions. idf and jdf are the (i,j) dimensions of ! domain on this processor, kmax is the number of model layers. !------------------------------------------------------------------- kmax = size(pref(:)) - 1 @@ -356,6 +381,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) idf = size(lonb) - 1 allocate( source_level(idf,jdf) ) + allocate( damp_level(idf,jdf) ) allocate( source_amp(idf,jdf) ) ! allocate( lat(idf,jdf) ) @@ -365,12 +391,16 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) ! ied as the source location via namelist input. !-------------------------------------------------------------------- do k=1,kmax + if (pref(k) < damp_level_pressure) then + klevel_of_damp = k + endif if (pref(k) > source_level_pressure) then klevel_of_source = k exit endif end do + do j=1,jdf !mj change of dimensions ! lat(:,j)= 0.5*( latb(:,j+1)+latb(:,j) ) @@ -378,12 +408,34 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) lat(i,j)= 0.5*( latb(j+1)+latb(j) ) source_level(i,j) = (kmax + 1) - ((kmax + 1 - & klevel_of_source)*cos(lat(i,j)) + 0.5) - source_amp(i,j) = Bt_0 + & - Bt_nh*0.5*(1.+tanh((lat(i,j)/pif-phi0n)/dphin)) + & - Bt_sh*0.5*(1.+tanh((lat(i,j)/pif-phi0s)/dphis)) + + damp_level(i,j) = klevel_of_damp !cig + +!code added by ipw - nov 23, 2016 + if (lat(i,j) > phi0n) then + source_amp(i,j) = Bt_0 + Bt_nh*0.5*(1.+tanh((lat(i,j)/pif-phi0n)/dphin))+ & + Bt_sh*0.5*(1.+tanh((lat(i,j)/pif-phi0s)/dphis)); + elseif (lat(i,j) < phi0s) then + source_amp(i,j) = Bt_0 + Bt_nh*0.5*(1.+tanh((lat(i,j)/pif-phi0n)/dphin))+ & + Bt_sh*0.5*(1.+tanh((lat(i,j)/pif-phi0s)/dphis)); + elseif ((lat(i,j) <= dphin) .and. (lat(i,j) >= dphis)) then + source_amp(i,j) = Bt_eq + elseif ((lat(i,j) <= phi0n) .and. (lat(i,j) > dphin)) then + source_amp(i,j) = Bt_0 + (Bt_eq-Bt_0)/(phi0n-dphin)*(phi0n-lat(i,j)) + elseif ((lat(i,j) < dphis) .and. (lat(i,j) >= phi0s)) then + source_amp(i,j) = Bt_0 + (Bt_eq-Bt_0)/(phi0s-dphis)*(phi0s-lat(i,j)) + endif + +! source_amp(i,j) = Bt_0 + & +! Bt_nh*0.5*(1.+tanh((lat(i,j)/pif-phi0n)/dphin)) + & +! Bt_sh*0.5*(1.+tanh((lat(i,j)/pif-phi0s)/dphis)) end do end do source_level = MIN (source_level, kmax-1) + damp_level = MIN (damp_level, kmax) + +!cig: make sure everyhing is ok +! write (*,*) "damp",pref(klevel_of_damp), ' ', klevel_of_damp, ' ', damp_level_pressure, ' ', damp_level(2,2), ' ', damp_level(12,2) ! deallocate( lat ) @@ -398,7 +450,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) !--------------------------------------------------------------------- ! if column diagnostics are desired, check that array dimensions are -! sufficiently large for the number of requests. +! sufficiently large for the number of requests. !--------------------------------------------------------------------- #ifdef COL_DIAG if (column_diagnostics_desired) then @@ -409,7 +461,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) endif !--------------------------------------------------------------------- -! allocate arrays needed for column diagnostics. +! allocate arrays needed for column diagnostics. !--------------------------------------------------------------------- allocate (do_column_diagnostics (jdf) ) allocate (diag_units (num_diag_pts) ) @@ -419,15 +471,15 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) allocate (diag_j (num_diag_pts) ) !--------------------------------------------------------------------- -! call initialize_diagnostic_columns to determine the locations -! (i, j, lat and lon) of any diagnostic columns in this processsor's +! call initialize_diagnostic_columns to determine the locations +! (i, j, lat and lon) of any diagnostic columns in this processsor's ! space and to open output files for the diagnostics. !--------------------------------------------------------------------- call initialize_diagnostic_columns & (mod_name, num_diag_pts_latlon, num_diag_pts_ij, & i_coords_gl, j_coords_gl, lat_coords_gl, & !mj dimensions are different -! lon_coords_gl, lonb(:,1), latb(1,:), do_column_diagnostics, & + lon_coords_gl, lonb(:,1), latb(1,:), do_column_diagnostics, & lon_coords_gl, lonb, latb, do_column_diagnostics, & diag_lon, diag_lat, diag_i, diag_j, diag_units) endif @@ -443,11 +495,11 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) do n=1,nc c0(n) = (n-1)*dc - cmax end do - + !-------------------------------------------------------------------- -! define the wavenumber kwv and its square k2 for the gravity waves -! contained in the spectrum. currently nk = 1, which means that the -! wavelength of all gravity waves considered is 300 km. +! define the wavenumber kwv and its square k2 for the gravity waves +! contained in the spectrum. currently nk = 1, which means that the +! wavelength of all gravity waves considered is 300 km. !-------------------------------------------------------------------- allocate ( kwv(nk) ) allocate ( k2 (nk) ) @@ -479,13 +531,20 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) register_diag_field (mod_name, 'kedy_cgwd', axes(1:3), Time, & 'effective eddy viscosity from cg_drag', 'm^2/s', & missing_value=missing_value) + id_source_level = & + register_diag_field (mod_name, 'source_level', axes(1:2), Time, & + 'effective eddy viscosity from cg_drag', 'm^2/s', & + missing_value=missing_value) !-------------------------------------------------------------------- -! allocate and define module variables to hold values across +! allocate and define module variables to hold values across ! timesteps, in the event that cg_drag is not called on every step. !-------------------------------------------------------------------- - allocate ( gwd_u(idf,jdf,kmax) ) - allocate ( gwd_v(idf,jdf,kmax) ) + +! allocate ( gwd_u(idf,jdf,kmax) ) RC 06/12/19 +! allocate ( gwd_v(idf,jdf,kmax) ) RC 06/12/19 + allocate ( gwd_u(is:ie,js:je,kmax) ) + allocate ( gwd_v(is:ie,js:je,kmax) ) !-------------------------------------------------------------------- ! if present, read the restart data file. @@ -516,12 +575,12 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) gwd_v(:,:,:) = 0.0 if (cg_drag_offset > 0) then cgdrag_alarm = cg_drag_offset - else + else cgdrag_alarm = cg_drag_freq endif !!$ endif !!$ vers = restart_versions(size(restart_versions(:))) -!!$ old_time_step = cgdrag_alarm +!!$ old_time_step = cgdrag_alarm !--------------------------------------------------------------------- ! mark the module as initialized. !--------------------------------------------------------------------- @@ -535,7 +594,7 @@ end subroutine cg_drag_init !#################################################################### - + subroutine cg_drag_time_vary (delt) real , intent(in) :: delt @@ -546,17 +605,17 @@ subroutine cg_drag_time_vary (delt) cgdrag_alarm = cgdrag_alarm - delt !--------------------------------------------------------------------- - + end subroutine cg_drag_time_vary !#################################################################### - + subroutine cg_drag_endts - + !-------------------------------------------------------------------- -! if this was a calculation step, reset cgdrag_alarm to indicate -! the time remaining before the next calculation of gravity wave +! if this was a calculation step, reset cgdrag_alarm to indicate +! the time remaining before the next calculation of gravity wave ! forcing. !-------------------------------------------------------------------- if (cgdrag_alarm <= 0 ) then @@ -570,9 +629,9 @@ end subroutine cg_drag_endts subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & Time, delt, gwfcng_x, gwfcng_y) -!-------------------------------------------------------------------- +!-------------------------------------------------------------------- ! cg_drag_calc defines the arrays needed to calculate the convective -! gravity wave forcing, calls gwfc to calculate the forcing, returns +! gravity wave forcing, calls gwfc to calculate the forcing, returns ! the desired output fields, and saves the values for later retrieval ! if they are not calculated on every timestep. ! @@ -589,7 +648,7 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & !------------------------------------------------------------------- ! intent(in) variables: ! -! is,js starting subdomain i,j indices of data in +! is,js starting subdomain i,j indices of data in ! the physics_window being integrated ! lat array of model latitudes at cell boundaries [radians] ! pfull pressure at model full levels [ Pa ] @@ -632,24 +691,25 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & ! local variables: ! ! dtdz temperature lapse rate [ deg K/m ] -! ked_gwfc effective diffusion coefficient from cg_drag_mod +! ked_gwfc effective diffusion coefficient from cg_drag_mod ! [ m^2/s ] ! zzchm heights at model levels [ m ] ! zu zonal velocity [ m/s ] ! zden atmospheric density [ kg/m^3 ] ! zbf buoyancy frequency [ /s ] -! gwd_xtnd zonal wind tendency resulting from cg_drag_mod +! gwd_xtnd zonal wind tendency resulting from cg_drag_mod ! [ m/s^2 ] -! ked_xtnd effective diffusion coefficient from cg_drag_mod +! ked_xtnd effective diffusion coefficient from cg_drag_mod ! [ m^2/s ] ! source_level k index of gravity wave source level ((i,j) array) +! damp_level k index of gravity wave mesospheric dumping level ((i,j) array) ! iz0 k index of gravity wave source level in a column ! used return code for netcdf diagnostics -! bflim minimum allowable value of squared buoyancy +! bflim minimum allowable value of squared buoyancy ! frequency [ /s^2 ] -! ie, je ending subdomain indices of data in the current +! ie, je ending subdomain indices of data in the current ! physics window being integrated -! imax, jmax, kmax +! imax, jmax, kmax ! physics window dimensions ! i, j, k, nn do loop indices ! @@ -665,24 +725,25 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & je = js + jmax - 1 !--------------------------------------------------------------------- -! if the convective gravity wave forcing should be calculated on +! if the convective gravity wave forcing should be calculated on ! this timestep (i.e., the alarm has gone off), proceed with the ! calculation. !--------------------------------------------------------------------- if (cgdrag_alarm <= 0) then - + write(6,*) 'GWD Calculated', cgdrag_alarm !----------------------------------------------------------------------- -! calculate temperature lapse rate. do one-sided differences over -! delta z at upper boundary and centered differences over 2 delta z +! calculate temperature lapse rate. do one-sided differences over +! delta z at upper boundary and centered differences over 2 delta z ! in the interior. dtdz is not needed at the lower boundary, since ! the source level is constrained to be above level kmax. !---------------------------------------------------------------------- do j=1,jmax do i=1,imax ! The following index-offsets are needed in case a physics_window is being used. - iz0 = source_level(i +is-1,j+js-1) + !iz0 = source_level(i +is-1,j+js-1) + iz0 = source_level(i,j) dtdz(i,j,1) = (temp (i,j,1) - temp (i,j,2))/ & (zfull(i,j,1) - zfull(i,j,2)) do k=2,iz0 @@ -698,14 +759,14 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & end do !---------------------------------------------------------------------- -! calculate buoyancy frequency. restrict the squared buoyancy +! calculate buoyancy frequency. restrict the squared buoyancy ! frequency to be no smaller than bflim. !---------------------------------------------------------------------- - do k=1,iz0 + do k=1,iz0 zbf(i,j,k) = (GRAV/temp(i,j,k))*(dtdz(i,j,k) + GRAV/CP_AIR) if (zbf(i,j,k) < bflim) then zbf(i,j,k) = sqrt(bflim) - else + else zbf(i,j,k) = sqrt(zbf(i,j,k)) endif end do @@ -725,7 +786,7 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & do k=1,iz0+1 zzchm(i,j,k) = zfull(i,j,k) end do - do k=1,iz0 + do k=1,iz0 zu(i,j,k) = uuu(i,j,k) zv(i,j,k) = vvv(i,j,k) end do @@ -734,8 +795,8 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & ! add an extra level above model top so that the gravity wave forcing ! occurring between the topmost model level and the upper boundary ! may be calculated. define variable values at the new top level as -! follows: z - use delta z of layer just below; u - extend vertical -! gradient occurring just below; density - geometric mean; buoyancy +! follows: z - use delta z of layer just below; u - extend vertical +! gradient occurring just below; density - geometric mean; buoyancy ! frequency - constant across model top. !---------------------------------------------------------------------- zzchm(i,j,0) = zzchm(i,j,1) + zzchm(i,j,1) - zzchm(i,j,2) @@ -745,25 +806,25 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & zbf(i,j,0) = zbf(i,j,1) end do end do - + !--------------------------------------------------------------------- ! pass the vertically-extended input arrays to gwfc. gwfc will cal- -! culate the gravity-wave forcing and, if desired, an effective eddy +! culate the gravity-wave forcing and, if desired, an effective eddy ! diffusion coefficient at each level above the source level. output ! is returned in the vertically-extended arrays gwfcng and ked_gwfc. -! upon return move the output fields into model-sized arrays. +! upon return move the output fields into model-sized arrays. !--------------------------------------------------------------------- - call gwfc (is, ie, js, je, source_level, source_amp, & + call gwfc (is, ie, js, je, damp_level, source_level, source_amp, lat, & zden, zu, zbf,zzchm, gwd_xtnd, ked_xtnd) gwfcng_x (:,:,1:kmax) = gwd_xtnd(:,:,1:kmax ) ked_gwfc_x(:,:,1:kmax) = ked_xtnd(:,:,1:kmax ) - - call gwfc (is, ie, js, je, source_level, source_amp, & + + call gwfc (is, ie, js, je, damp_level, source_level, source_amp, lat, & zden, zv, zbf,zzchm, gwd_ytnd, ked_ytnd) gwfcng_y (:,:,1:kmax) = gwd_ytnd(:,:,1:kmax ) ked_gwfc_y(:,:,1:kmax) = ked_ytnd(:,:,1:kmax ) - + !-------------------------------------------------------------------- ! store the gravity wave forcing into a processor-global array. !------------------------------------------------------------------- @@ -771,91 +832,89 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & gwd_v(is:ie,js:je,:) = gwfcng_y(:,:,:) -#ifdef COL_DIAG -!-------------------------------------------------------------------- -! if column diagnostics are desired, determine if any columns are on -! this processor. if so, call column_diagnostics_header to write -! out location and timestamp information. then output desired -! quantities to the diag_unit file. -!--------------------------------------------------------------------- - if (column_diagnostics_desired) then - do j=1,jmax - if (do_column_diagnostics(j+js-1)) then - do nn=1,num_diag_pts - if (js + j - 1 == diag_j(nn)) then - call column_diagnostics_header & - (mod_name, diag_units(nn), Time, nn, diag_lon, & - diag_lat, diag_i, diag_j) - iz0 = source_level (diag_i(nn), j) - write (diag_units(nn),'(a, i5)') & - ' source_level =', iz0 - write (diag_units(nn),'(a)') & - ' k u z density& - & bf gwforcing' - do k=0,iz0 - write (diag_units(nn), '(i5, 2x, 5e12.5)') & - k, & - zu (diag_i(nn),j,k), & - zzchm (diag_i(nn),j,k), & - zden (diag_i(nn),j,k), & - zbf (diag_i(nn),j,k), & - gwd_xtnd (diag_i(nn),j,k) - end do - write (diag_units(nn), '(i5, 14x, 2e12.5)') & - iz0+1, & - zzchm (diag_i(nn),j,iz0+1), & - zden (diag_i(nn),j,iz0+1) - endif - end do ! (nn loop) - endif ! (do_column_diagnostics) - end do ! (j loop) - endif ! (column_diagnostics_desired) -#endif +! #ifdef COL_DIAG +! !-------------------------------------------------------------------- +! ! if column diagnostics are desired, determine if any columns are on +! ! this processor. if so, call column_diagnostics_header to write +! ! out location and timestamp information. then output desired +! ! quantities to the diag_unit file. +! !--------------------------------------------------------------------- +! if (column_diagnostics_desired) then +! do j=1,jmax +! if (do_column_diagnostics(j+js-1)) then +! do nn=1,num_diag_pts +! if (js + j - 1 == diag_j(nn)) then +! call column_diagnostics_header & +! (mod_name, diag_units(nn), Time, nn, diag_lon, & +! diag_lat, diag_i, diag_j) +! iz0 = source_level (diag_i(nn), j) +! write (diag_units(nn),'(a, i5)') & +! ' source_level =', iz0 +! write (diag_units(nn),'(a)') & +! ' k u z density& +! & bf gwforcing' +! do k=0,iz0 +! write (diag_units(nn), '(i5, 2x, 5e12.5)') & +! k, & +! zu (diag_i(nn),j,k), & +! zzchm (diag_i(nn),j,k), & +! zden (diag_i(nn),j,k), & +! zbf (diag_i(nn),j,k), & +! gwd_xtnd (diag_i(nn),j,k) +! end do +! write (diag_units(nn), '(i5, 14x, 2e12.5)') & +! iz0+1, & +! zzchm (diag_i(nn),j,iz0+1), & +! zden (diag_i(nn),j,iz0+1) +! endif +! end do ! (nn loop) +! endif ! (do_column_diagnostics) +! end do ! (j loop) +! endif ! (column_diagnostics_desired) +! #endif !-------------------------------------------------------------------- -! if activated, store the effective eddy diffusivity into a -! processor-global array, and if desired as a netcdf diagnostic, +! if activated, store the effective eddy diffusivity into a +! processor-global array, and if desired as a netcdf diagnostic, ! send the data to diag_manager_mod. !------------------------------------------------------------------- if (id_kedx_cgwd > 0) then -! used = send_data (id_kedx_cgwd, ked_gwfc_x, Time, is, js, 1) - used = send_data (id_kedx_cgwd, ked_gwfc_x, Time) !st 2013 FMS seems to have paralelisation issues when called with ...Time, is, js, 1). + used = send_data (id_kedx_cgwd, ked_gwfc_x, Time)!, is, js, 1) endif if (id_kedy_cgwd > 0) then -! used = send_data (id_kedy_cgwd, ked_gwfc_y, Time, is, js, 1) - used = send_data (id_kedy_cgwd, ked_gwfc_y, Time) + used = send_data (id_kedy_cgwd, ked_gwfc_y, Time)!, is, js, 1) endif - + if (id_source_level > 0) then + used = send_data (id_source_level, real(source_level), Time) + endif !-------------------------------------------------------------------- ! save any other netcdf file diagnostics that are desired. !-------------------------------------------------------------------- if (id_bf_cgwd > 0) then -! used = send_data (id_bf_cgwd, zbf(:,:,1:), Time, is, js ) - used = send_data (id_bf_cgwd, zbf(:,:,1:), Time) + used = send_data (id_bf_cgwd, zbf(:,:,1:), Time)!, is, js ) endif if (id_gwfx_cgwd > 0) then -! used = send_data (id_gwfx_cgwd, gwfcng_x, Time, is, js, 1) - used = send_data (id_gwfx_cgwd, gwfcng_x, Time) + used = send_data (id_gwfx_cgwd, gwfcng_x, Time)!, is, js, 1) endif if (id_gwfy_cgwd > 0) then -! used = send_data (id_gwfy_cgwd, gwfcng_y, Time, is, js, 1) - used = send_data (id_gwfy_cgwd, gwfcng_y, Time) + used = send_data (id_gwfy_cgwd, gwfcng_y, Time)!, is, js, 1) endif !-------------------------------------------------------------------- -! if this is not a timestep on which gravity wave forcing is to be +! if this is not a timestep on which gravity wave forcing is to be ! calculated, retrieve the values calculated previously from storage ! and return to the calling subroutine. !-------------------------------------------------------------------- else ! (cgdrag_alarm <= 0) + write(6,*) 'GWD not recalculated',cgdrag_alarm gwfcng_x(:,:,:) = gwd_u(is:ie,js:je,:) gwfcng_y(:,:,:) = gwd_v(is:ie,js:je,:) endif ! (cgdrag_alarm <= 0) @@ -865,7 +924,7 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & ! will recalculate the NOGWD tendencies call cg_drag_endts call cg_drag_time_vary(delt) - + end subroutine cg_drag_calc @@ -927,8 +986,8 @@ end subroutine cg_drag_end !!$ unit = open_restart_file ('RESTART/cg_drag.res', 'write') !!$ !!$!------------------------------------------------------------------- -!!$! the root pe writes out the restart version, the time remaining -!!$! before the next call to cg_drag_mod and the current cg_drag +!!$! the root pe writes out the restart version, the time remaining +!!$! before the next call to cg_drag_mod and the current cg_drag !!$! timestep. !!$!------------------------------------------------------------------- !!$ if (mpp_pe() == mpp_root_pe() ) then @@ -937,7 +996,7 @@ end subroutine cg_drag_end !!$ endif !!$ !!$!------------------------------------------------------------------- -!!$! each processor writes out its gravity wave forcing tendency +!!$! each processor writes out its gravity wave forcing tendency !!$! on the zonal flow. !!$!------------------------------------------------------------------- !!$ call write_data (unit, gwd_u) @@ -968,11 +1027,11 @@ end subroutine cg_drag_end !!$ real :: secs_per_day = SECONDS_PER_DAY !!$ !!$!------------------------------------------------------------------- -!!$! local variables: -!!$! -!!$! unit unit number for nml file -!!$! chvers character representation of restart version -!!$! vers restart version +!!$! local variables: +!!$! +!!$! unit unit number for nml file +!!$! chvers character representation of restart version +!!$! vers restart version !!$! dummy array to hold restart version 1 control variables !!$! old_time_step cg_drag timestep used in previous model run [ s ] !!$! secs_per_day seconds in a day [ s ] @@ -981,7 +1040,7 @@ end subroutine cg_drag_end !!$ !!$ !!$!-------------------------------------------------------------------- -!!$! open file to read restart data. +!!$! open file to read restart data. !!$!--------------------------------------------------------------------- !!$ unit = open_restart_file ('INPUT/cg_drag.res','read') !!$ @@ -997,36 +1056,36 @@ end subroutine cg_drag_end !!$ endif !!$ !!$!-------------------------------------------------------------------- -!!$! read control information from restart file. +!!$! read control information from restart file. !!$!-------------------------------------------------------------------- !!$ if (vers == 1) then !!$ !!$!-------------------------------------------------------------------- !!$! if reading restart version 1, use the contents of array dummy to -!!$! define the cg_drag timestep that was used in the run which wrote -!!$! the restart. define the time remaining before the next cg_drag +!!$! define the cg_drag timestep that was used in the run which wrote +!!$! the restart. define the time remaining before the next cg_drag !!$! calculation to either be the previous timestep or the current !!$! offset, if that is specified. this assumes that the restart was !!$! written at 00Z. !!$!-------------------------------------------------------------------- -!!$ read (unit) dummy +!!$ read (unit) dummy !!$ old_time_step = secs_per_day*dummy(4) + dummy(3) !!$ if (cg_drag_offset == 0) then !!$ cgdrag_alarm = old_time_step !!$ else -!!$ cgdrag_alarm = cg_drag_offset +!!$ cgdrag_alarm = cg_drag_offset !!$ endif -!!$ else +!!$ else !!$ !!$!-------------------------------------------------------------------- -!!$! for restart version 2, read the time remaining until the next +!!$! for restart version 2, read the time remaining until the next !!$! cg_drag calculation, and the previously used timestep. !!$!--------------------------------------------------------------------- !!$ read (unit) cgdrag_alarm, old_time_step !!$ endif !!$ !!$!------------------------------------------------------------------- -!!$! read restart data (gravity wave forcing tendency terms) and close +!!$! read restart data (gravity wave forcing tendency terms) and close !!$! unit. !!$!------------------------------------------------------------------- !!$ call read_data (unit, gwd_u) @@ -1034,8 +1093,8 @@ end subroutine cg_drag_end !!$ call close_file (unit) !!$ !!$!-------------------------------------------------------------------- -!!$! if current cg_drag calling frequency differs from that previously -!!$! used, adjust the time remaining before the next calculation. +!!$! if current cg_drag calling frequency differs from that previously +!!$! used, adjust the time remaining before the next calculation. !!$!-------------------------------------------------------------------- !!$ if (cg_drag_freq /= old_time_step) then !!$ cgdrag_alarm = cgdrag_alarm - old_time_step + cg_drag_freq @@ -1048,7 +1107,7 @@ end subroutine cg_drag_end !!$ !!$!-------------------------------------------------------------------- !!$! if cg_drag_offset is specified and is smaller than the time remain- -!!$! ing until the next calculation, modify the time remaining to be +!!$! ing until the next calculation, modify the time remaining to be !!$! that offset time. the assumption is made that the restart was !!$! written at 00Z. !!$!-------------------------------------------------------------------- @@ -1066,8 +1125,8 @@ end subroutine cg_drag_end !!$ !!$subroutine read_nc_restart_file !!$!----------------------------------------------------------------------- -!!$! subroutine read_restart_nc reads a netcdf restart file to obtain -!!$! the variables needed upon experiment restart. +!!$! subroutine read_restart_nc reads a netcdf restart file to obtain +!!$! the variables needed upon experiment restart. !!$!----------------------------------------------------------------------- !!$ !!$!--------------------------------------------------------------------- @@ -1111,8 +1170,8 @@ end subroutine cg_drag_end !!$ vers = restart_versions(size(restart_versions(:))) !!$ !!$!-------------------------------------------------------------------- -!!$! if current cg_drag calling frequency differs from that previously -!!$! used, adjust the time remaining before the next calculation. +!!$! if current cg_drag calling frequency differs from that previously +!!$! used, adjust the time remaining before the next calculation. !!$!-------------------------------------------------------------------- !!$ if (cg_drag_freq /= old_time_step) then !!$ cgdrag_alarm = cgdrag_alarm - old_time_step + cg_drag_freq @@ -1126,7 +1185,7 @@ end subroutine cg_drag_end !!$ !!$!-------------------------------------------------------------------- !!$! if cg_drag_offset is specified and is smaller than the time remain- -!!$! ing until the next calculation, modify the time remaining to be +!!$! ing until the next calculation, modify the time remaining to be !!$! that offset time. the assumption is made that the restart was !!$! written at 00Z. !!$!-------------------------------------------------------------------- @@ -1144,10 +1203,10 @@ end subroutine cg_drag_end !!$subroutine cg_drag_register_restart !!$ !!$ character(len=64) :: fname = 'cg_drag.res.nc' ! name of restart file -!!$ character(len=64) :: fname2 +!!$ character(len=64) :: fname2 !!$ integer :: id_restart !!$ -!!$ call get_mosaic_tile_file(fname, fname2, .false. ) +!!$ call get_mosaic_tile_file(fname, fname2, .false. ) !!$ allocate(Cg_restart) !!$ if(trim(fname2) == trim(fname)) then !!$ Til_restart => Cg_restart @@ -1172,10 +1231,10 @@ end subroutine cg_drag_end !!$! !!$! !!$! write out restart file. -!!$! Arguments: -!!$! timestamp (optional, intent(in)) : A character string that represents the model time, +!!$! Arguments: +!!$! timestamp (optional, intent(in)) : A character string that represents the model time, !!$! used for writing restart. timestamp will append to -!!$! the any restart file name as a prefix. +!!$! the any restart file name as a prefix. !!$! !!$! !!$subroutine cg_drag_restart(timestamp) @@ -1190,21 +1249,21 @@ end subroutine cg_drag_end !#################################################################### -subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & +subroutine gwfc (is, ie, js, je, damp_level, source_level, source_amp, lat, rho, u, & bf, z, gwf, ked) !------------------------------------------------------------------- ! subroutine gwfc computes the gravity wave-driven-forcing on the -! zonal wind given vertical profiles of wind, density, and buoyancy -! frequency. -! Based on version implemented in SKYHI -- 27 Oct 1998 by M.J. +! zonal wind given vertical profiles of wind, density, and buoyancy +! frequency. +! Based on version implemented in SKYHI -- 27 Oct 1998 by M.J. ! Alexander and L. Bruhwiler. !------------------------------------------------------------------- !------------------------------------------------------------------- integer, intent(in) :: is, ie, js, je -integer, dimension(:,:), intent(in) :: source_level -real, dimension(:,:), intent(in) :: source_amp +integer, dimension(:,:), intent(in) :: source_level, damp_level +real, dimension(:,:), intent(in) :: source_amp, lat real, dimension(:,:,0:), intent(in) :: rho, u, bf, z real, dimension(:,:,0:), intent(out) :: gwf real, dimension(:,:,0:), intent(out) :: ked @@ -1216,9 +1275,10 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ! in the physics_window being integrated ! source_level k index of model level serving as gravity wave ! source +! damp_level k index of the lowest model level at which all drag that reaches the model top is partially dumped ! source_amp amplitude of gravity wave source [Pa] -! -! rho atmospheric density [ kg/m^3 ] +! +! rho atmospheric density [ kg/m^3 ] ! u zonal wind component [ m/s ] ! bf buoyancy frequency [ /s ] ! z height of model levels [ m ] @@ -1229,7 +1289,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ! ! intent(out), optional variables: ! -! ked eddy diffusion coefficient from gravity wave +! ked eddy diffusion coefficient from gravity wave ! forcing [ m^2/s ] ! !------------------------------------------------------------------ @@ -1244,28 +1304,28 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & real , dimension (nc) :: c0mu0, B0 real :: fm, fe, Hb, alp2, Foc, c, test, rbh,& eps, Bsum - integer :: iz0 + integer :: iz0, iztop integer :: i, j, k, ink, n - real :: ampl + real :: ampl, cwthis, Bnthis, flagthis !------------------------------------------------------------------ ! local variables: -! +! ! wv_frcng gravity wave forcing tendency [ m/s^2 ] ! diff_coeff eddy diffusion coefficient [ m2/s ] -! c0mu difference between phase speed of wave n and u +! c0mu difference between phase speed of wave n and u ! [ m/s ] ! dz delta z between model levels [ m ] -! fac factor used in determining if wave is breaking +! fac factor used in determining if wave is breaking ! [ s/m ] -! omc critical frequency that marks total internal +! omc critical frequency that marks total internal ! reflection [ /s ] -! msk indicator as to whether wave n is still propagating -! upwards (msk=1), or has been removed from the +! msk indicator as to whether wave n is still propagating +! upwards (msk=1), or has been removed from the ! spectrum because of breaking or reflection (msk=0) ! c0mu0 difference between phase speed of wave n and u at the ! source level [ m/s ] ! B0 wave momentum flux amplitude for wave n [ (m/s)^2 ] -! fm used to sum up momentum flux from all waves n +! fm used to sum up momentum flux from all waves n ! deposited at a level [ (m/s)^2 ] ! fe used to sum up contributions to diffusion coefficient ! from all waves n at a level [ (m/s)^3 ] @@ -1278,18 +1338,18 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ! rbh atmospheric density at half-level (geometric mean) ! [ kg/m^3 ] ! eps intermittency factor -! Bsum total mag of gravity wave momentum flux at source +! Bsum total mag of gravity wave momentum flux at source ! level, divided by the density [ m^2/s^2 ] ! iz0 source level vertical index for the given column ! i,j,k spatial do loop indices -! ink wavenumber loop index -! n phase speed loop index -! ampl gravity wave stress [Pa] +! ink wavenumber loop index +! n phase speed loop index +! ampl gravity wave stress [Pa] ! !-------------------------------------------------------------------- !------------------------------------------------------------------- -! initialize the output arrays. these will hold values at each +! initialize the output arrays. these will hold values at each ! (i,j,k) point, summed over the wavelengths and phase speeds ! defining the gravity wave spectrum. !------------------------------------------------------------------- @@ -1297,19 +1357,39 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ked = 0.0 do j=1,size(u,2) - do i=1,size(u,1) -! The following index-offsets are needed in case a physics_window is being used. - iz0 = source_level(i+is-1,j+js-1) - ampl= source_amp(i+is-1,j+js-1) + + do i=1,size(u,1) +!added by cig, january 2017, Changed indicies RC 2019 + !if ((lat(i+is-1,j+js-1) <= dphin) .and. (lat(i+is-1,j+js-1) >= dphis)) then + if ((lat(i,j) <= dphin) .and. (lat(i,j) >= dphis)) then + cwthis=cwtropics + Bnthis=0. + flag=0 + else + cwthis=cw + Bnthis=Bn + flagthis=flag + endif + +! RC Changed Indicies The following index-offsets are needed in case a physics_window is being used. +! iz0 = source_level(i+is-1,j+js-1) +! iztop = damp_level(i+is-1,j+js-1) +! ampl= source_amp(i+is-1,j+js-1) + +! The following index-offsets are needed in case a physics_window is being used. + iz0 = source_level(i,j) + ampl= source_amp(i,j) +! DO I ADD iztop=damp_level(i,j)? + iztop=damp_level(i,j) !-------------------------------------------------------------------- -! define wave momentum flux (B0) at source level for each phase -! speed n, and the sum over all phase speeds (Bsum), which is needed -! to calculate the intermittency. +! define wave momentum flux (B0) at source level for each phase +! speed n, and the sum over all phase speeds (Bsum), which is needed +! to calculate the intermittency. !------------------------------------------------------------------- Bsum = 0. do n=1,nc - c0mu0(n) = c0(n) - u(i,j,iz0) + c0mu0(n) = c0(n) - u(i,j,iz0) !--------------------------------------------------------------------- ! when the wave phase speed is same as wind speed, there is no @@ -1317,19 +1397,21 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & !--------------------------------------------------------------------- if (c0mu0(n) == 0.0) then B0(n) = 0.0 - else + else !--------------------------------------------------------------------- ! define wave momentum flux at source level for phase speed n. Add ! the contribution from this phase speed to the previous sum. !--------------------------------------------------------------------- - c = c0(n)*flag + c0mu0(n)*(1 - flag) + c = c0(n)*flagthis + c0mu0(n)*(1 - flagthis) if (c0mu0(n) < 0.0) then - B0(n) = -1.0*(Bw*exp(-alog(2.0)*(c/cw)**2) + & - Bn*exp(-alog(2.0)*(c/cn)**2)) - else - B0(n) = (Bw*exp(-alog(2.0)*(c/cw)**2) + & - Bn*exp(-alog(2.0)*(c/cn)**2)) + B0(n) = -1.0*(Bw*exp(-alog(2.0)*(c/cwthis)**2) + & + Bnthis*exp(-alog(2.0)*(c/cn)**2)) + B0(n) =B0(n)*kelvin_kludge + else + B0(n) = (Bw*exp(-alog(2.0)*(c/cwthis)**2) + & + Bnthis*exp(-alog(2.0)*(c/cn)**2)) + endif Bsum = Bsum + abs(B0(n)) endif @@ -1364,7 +1446,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & fac(k) = 0.5*(rho(i,j,k)/rho(i,j,iz0))*kwv(ink)/bf(i,j,k) end do - do k=0,iz0 + do k=0,iz0 dz(k) = z(i,j,k) - z(i,j,k+1) Hb = -(dz(k))/alog(rho(i,j,k)/rho(i,j,k+1)) alp2 = 0.25/(Hb*Hb) @@ -1373,14 +1455,14 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & end do !--------------------------------------------------------------------- -! initialize a flag which will indicate which waves are still +! initialize a flag which will indicate which waves are still ! propagating upwards. !--------------------------------------------------------------------- msk = 1 !---------------------------------------------------------------------- -! integrate upwards from the source level. define variables over -! which to sum the deposited flux and effective eddy diffusivity +! integrate upwards from the source level. define variables over +! which to sum the deposited flux and effective eddy diffusivity ! from all waves breaking at a given level. !---------------------------------------------------------------------- do k=iz0, 0, -1 @@ -1395,7 +1477,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & c0mu(k) = c0(n) - u(i,j,k) !---------------------------------------------------------------------- -! if phase speed matches the wind speed, remove c0(n) from the +! if phase speed matches the wind speed, remove c0(n) from the ! set of propagating waves. ! epg: This seems to be an unphysical decision, as the wave should ! break, having reached a critical level. But it's extremely @@ -1408,7 +1490,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & else !--------------------------------------------------------------------- -! define the criterion which determines if wave is reflected at this +! define the criterion which determines if wave is reflected at this ! level (test). !--------------------------------------------------------------------- test = abs(c0mu(k))*kwv(ink) - omc(k) @@ -1419,22 +1501,22 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ! propagating set. !--------------------------------------------------------------------- msk(n) = 0 - else + else !--------------------------------------------------------------------- -! if wave is not reflected at this level, determine if it is -! breaking at this level (Foc >= 0), or if wave speed relative to -! windspeed has changed sign from its value at the source level +! if wave is not reflected at this level, determine if it is +! breaking at this level (Foc >= 0), or if wave speed relative to +! windspeed has changed sign from its value at the source level ! (c0mu0(n)*c0mu <= 0). if it is above the source level and is -! breaking, then add its momentum flux to the accumulated sum at -! this level, and increase the effective diffusivity accordingly. +! breaking, then add its momentum flux to the accumulated sum at +! this level, and increase the effective diffusivity accordingly. ! set flag to remove phase speed c0(n) from the set of active waves ! moving upwards to the next level. !--------------------------------------------------------------------- -! epg: if you are at the model top, deposit all momentum here, to +! epg: if you are at the model top, deposit all momentum here, to ! prevent waves from escaping the top of the model. See -! Shaw et al. 2010? for details on why this in important. +! Shaw et al. 2010? for details on why this in important. if ( k==0 ) then msk(n) = 0 if (k < iz0) then @@ -1442,7 +1524,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & fe = fe + c0mu(k)*B0(n) endif else - + Foc = B0(n)/(c0mu(k) )**3 - fac(k) if ((Foc >= 0.0) .or. & (c0mu0(n)*c0mu(k) <= 0.0)) then @@ -1459,51 +1541,77 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & end do ! phase speed loop !---------------------------------------------------------------------- -! compute the gravity wave momentum flux forcing and eddy +! compute the gravity wave momentum flux forcing and eddy ! diffusion coefficient obtained across the entire wave spectrum ! at this level. !---------------------------------------------------------------------- if ( k < iz0) then rbh = sqrt(rho(i,j,k)*rho(i,j,k+1)) wv_frcng(k) = ( rho(i,j,iz0)/rbh)*fm*eps/dz(k) - - !epg: enforce momentum conservation at model top; in this case, all the momentum + + !epg: enforce momentum conservation at model top; in this case, all the momentum ! deposited in the uppermost layer (which exist above the top model level, ! as explained in cg_drag_calc, must be added to the level below, which is ! the actual top level of the model. + !cig: place the extra momentum flux in the top 3 layers instead of all in the top layer if ( k==0 ) then - wv_frcng(k+1) = 0.5*wv_frcng(k+1) + wv_frcng(k) + wv_frcng(k+1) = 0.5*wv_frcng(k+1) !+ weighttop*wv_frcng(k) cig commented out else - wv_frcng(k+1) = 0.5*(wv_frcng(k+1) + wv_frcng(k)) + wv_frcng(k+1) = 0.5*(wv_frcng(k+1) + wv_frcng(k)) endif - + + + + diff_coeff(k) = (rho(i,j,iz0)/rbh)*fe*eps/(dz(k)* & bf(i,j,k)*bf(i,j,k)) - + !epg: following what we did above... + !cig: place the extra momentum flux in the top 3 layers instead of all in the top layer if ( k==0 ) then - diff_coeff(k+1) = 0.5*diff_coeff(k+1) + diff_coeff(k) + diff_coeff(k+1) = 0.5*diff_coeff(k+1) !+ weighttop*diff_coeff(k) cig commented out else diff_coeff(k+1) = 0.5*(diff_coeff(k+1) + diff_coeff(k)) endif - - else + + !cig: following what we did above... + + + else wv_frcng(iz0) = 0.0 diff_coeff(iz0) = 0.0 endif - end do ! (k loop) + end do ! (k loop) + +!cig: place the extra momentum flux in the layers above a specific threshold instead of all in the top layer (k=0 isn't a real model level) +! write (*,*) "iztop",iztop, ' ', damp_level(i,j), ' ', damp_level_pressure, ' ', iz0, ' ', source_level_pressure + do k=1,iztop + wv_frcng(k) = wv_frcng(k) + wv_frcng(0)/REAL(iztop) + diff_coeff(k) = diff_coeff(k) + diff_coeff(0)/REAL(iztop) + end do + +!cig: place the extra momentum flux in the top 3 layers instead of all in the top layer +! wv_frcng(1) = wv_frcng(1) + weighttop*wv_frcng(0) +! wv_frcng(2) = wv_frcng(2) + weightminus1*wv_frcng(0) +! wv_frcng(3) = wv_frcng(3) + weightminus2*wv_frcng(0) + +! diff_coeff(1) = diff_coeff(1) + weighttop*diff_coeff(0) +! diff_coeff(2) = diff_coeff(2) + weightminus1*diff_coeff(0) +! diff_coeff(3) = diff_coeff(3) + weightminus2*diff_coeff(0) !--------------------------------------------------------------------- ! increment the total forcing at each point with that obtained from ! the set of waves with the current wavenumber. !--------------------------------------------------------------------- - do k=0,iz0 + + + do k=0,iz0 gwf(i,j,k) = gwf(i,j,k) + wv_frcng(k) ked(i,j,k) = ked(i,j,k) + diff_coeff(k) - end do + end do end do ! wavelength loop - end do ! i loop - end do ! j loop + end do ! i loop + end do ! j loop !-------------------------------------------------------------------- @@ -1517,5 +1625,3 @@ end subroutine gwfc end module cg_drag_mod - - diff --git a/src/atmos_param/damping_driver/damping_driver.f90 b/src/atmos_param/damping_driver/damping_driver.f90 index f1829d5d3..6c3881595 100644 --- a/src/atmos_param/damping_driver/damping_driver.f90 +++ b/src/atmos_param/damping_driver/damping_driver.f90 @@ -44,7 +44,7 @@ module damping_driver_mod integer :: nlev_rayfric real :: sponge_pbottom = 50. ! [Pa] logical :: do_mg_drag = .false. -!epg: Use cg_drag.f90, GFDL's version of the Alexander and Dunkerton 1999 +!epg: Use cg_drag.f90, GFDL's version of the Alexander and Dunkerton 1999 ! Non-orographic gravity wave parameterization, updated as for Cohen et al. 2013 ! mj actively choose rayleigh friction logical :: do_rayleigh = .false. @@ -101,7 +101,7 @@ module damping_driver_mod real :: rfactr -! note: +! note: ! rfactr = coeff. for damping momentum at the top level character(len=128) :: version = '$Id: damping_driver.f90,v 10.0 2003/10/24 22:00:25 fms Exp $' @@ -119,7 +119,7 @@ subroutine damping_driver (is, js, lat, Time, delt, pfull, phalf, zfull, zhalf, u, v, t, q, r, udt, vdt, tdt, qdt, rdt, & ! mask, kbot) z_pbl, mask, kbot) - + !----------------------------------------------------------------------- integer, intent(in) :: is, js real, dimension(:,:), intent(in) :: lat @@ -265,6 +265,7 @@ subroutine damping_driver (is, js, lat, Time, delt, pfull, phalf, zfull, zhalf, !mj updating call to riga version of cg_drag !call cg_drag_calc (is, js, lat, pfull, zfull, t, u, Time, & ! delt, utnd) + call cg_drag_calc (is, js, lat, pfull, zfull, t, u, v, Time, delt, utnd, vtnd) udt = udt + utnd vdt = vdt + vtnd !mj @@ -274,9 +275,9 @@ subroutine damping_driver (is, js, lat, Time, delt, pfull, phalf, zfull, zhalf, if ( id_udt_cgwd > 0 ) then ! used = send_data ( id_udt_cgwd, utnd, Time, is, js, 1, & ! rmask=mask ) - used = send_data ( id_udt_cgwd, utnd, Time, rmask=mask ) + used = send_data ( id_udt_cgwd, utnd, Time) endif - + endif ! constant drag, modeled on Alexander-Dunkerton winter average @@ -298,7 +299,7 @@ subroutine damping_driver (is, js, lat, Time, delt, pfull, phalf, zfull, zhalf, *( -1.65*abs(lat)**3 +2.5*lat**2 +0.17*abs(lat) +const_drag_off ) end where enddo - udt = udt + utnd + udt = udt + utnd !----- diagnostics ----- @@ -366,12 +367,13 @@ end subroutine damping_driver !####################################################################### - subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) + subroutine damping_driver_init (is, ie, js, je, lonb, latb, pref, axes, Time, sgsmtn) real, intent(in) :: lonb(:), latb(:), pref(:) - integer, intent(in) :: axes(4) + integer, dimension(4), intent(in) :: axes type(time_type), intent(in) :: Time real, dimension(:,:), intent(out) :: sgsmtn + integer, intent(in) :: is, ie, js, je !----------------------------------------------------------------------- ! lonb = longitude in radians of the grid box edges ! latb = latitude in radians of the grid box edges @@ -426,13 +428,13 @@ subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) !----------------------------------------------------------------------- !----- mountain gravity wave drag ----- - if (do_mg_drag) call mg_drag_init (lonb, latb, sgsmtn) + if (do_mg_drag) call mg_drag_init (is, ie, js, je, lonb, latb, sgsmtn) !-------------------------------------------------------------------- !----- Alexander-Dunkerton gravity wave drag ----- - + if (do_cg_drag) then - call cg_drag_init (lonb, latb, pref, Time=Time, axes=axes) + call cg_drag_init (is, ie, js, je, lonb, latb, pref, Time=Time, axes=axes) endif !----------------------------------------------------------------------- @@ -454,7 +456,7 @@ subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) register_diag_field ( mod_name, 'tdt_diss_rdamp', axes(1:3), Time, & 'Dissipative heating from Rayleigh damping',& 'deg_k/s', missing_value=missing_value ) - + id_diss_heat_rdamp = & register_diag_field ( mod_name, 'diss_heat_rdamp', axes(1:2), Time, & 'Integrated dissipative heating from Rayleigh damping',& @@ -499,7 +501,7 @@ subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) register_diag_field ( mod_name, 'tdt_diss_gwd', axes(1:3), Time, & 'Dissipative heating from gravity wave drag',& 'deg_k/s', missing_value=missing_value ) - + id_diss_heat_gwd = & register_diag_field ( mod_name, 'diss_heat_gwd', axes(1:2), Time, & 'Integrated dissipative heating from gravity wave drag',& @@ -521,7 +523,7 @@ subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) 'u wind tendency for constant drag', 'm/s2', & missing_value=missing_value ) endif - + !----------------------------------------------------------------------- !----- topo wave drag ----- @@ -639,4 +641,3 @@ end subroutine rayleigh !####################################################################### end module damping_driver_mod - diff --git a/src/atmos_param/mg_drag/mg_drag.f90 b/src/atmos_param/mg_drag/mg_drag.f90 index 1279a852c..1e0ac7701 100644 --- a/src/atmos_param/mg_drag/mg_drag.f90 +++ b/src/atmos_param/mg_drag/mg_drag.f90 @@ -6,7 +6,7 @@ module mg_drag_mod !------------------------------------------------------------------- ! Calculates partial tendencies for the zonal and meridional winds -! due to the effect of mountain gravity wave drag +! due to the effect of mountain gravity wave drag !------------------------------------------------------------------- use topography_mod, only: get_topog_stdev @@ -47,13 +47,13 @@ module mg_drag_mod !--------------------------------------------------------------------- ! xl_mtn effective mountain length ( set currently to 100km) ! acoef order unity "tunable" parameter -! gmax order unity "tunable" parameter +! gmax order unity "tunable" parameter ! (may be enhanced to increase drag) ! rho stand value for density of the air at sea-level (1.13 KG/M**3) ! low_lev_frac - fraction of atmosphere (from bottom up) considered ! to be "low-level-layer for base flux calc. and where no ! wave breaking is allowed. -! flux_cut_level pressure level (Pa) above which flux divergence is set to zero +! flux_cut_level pressure level (Pa) above which flux divergence is set to zero !----------------------------------------------------------------------- real :: & @@ -65,22 +65,34 @@ module mg_drag_mod ,low_lev_frac = .23 real :: flux_cut_level= 0.0 +real :: tp_level = .9 + +integer :: ideal_lat_cent_NH = 48 +integer :: ideal_lat_cent_SH = 17 +integer :: ideal_pfull_cent = 19 +real :: ideal_strength = -0.00002 logical :: do_netcdf_restart = .true. logical :: do_conserve_energy = .false. logical :: do_mcm_mg_drag = .false. +logical :: no_strat_drag = .false. +logical :: no_trop_drag = .false. +logical :: do_idealised_NH = .false. +logical :: do_idealised_SH = .false. character(len=128) :: source_of_sgsmtn = 'input' namelist / mg_drag_nml / do_netcdf_restart, & xl_mtn, gmax, acoef, rho, low_lev_frac, & do_conserve_energy, do_mcm_mg_drag, & - source_of_sgsmtn, flux_cut_level + source_of_sgsmtn, flux_cut_level, & + no_strat_drag, no_trop_drag, tp_level, & + do_idealised_NH, do_idealised_SH, ideal_lat_cent_NH, ideal_lat_cent_SH, ideal_pfull_cent, ideal_strength public mg_drag, mg_drag_init, mg_drag_end contains -!############################################################################# +!############################################################################# subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & zfull,zhalf,dtaux,dtauy,dtemp,taubx, tauby, tausf,& @@ -129,14 +141,14 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & ! (dimensioned IDIM x JDIM)-kg/m/s**2 ! = -(RHO*U**3/(N*XL))*G(FR) FOR N**2 > 0 ! = 0 FOR N**2 <=0 -! DTAUX Tendency of the zonal wind component deceleration +! DTAUX Tendency of the zonal wind component deceleration ! (dimensioned IDIM x JDIM x KDIM) -! DTAUY Tendency of the meridional wind component deceleration +! DTAUY Tendency of the meridional wind component deceleration ! (dimensioned IDIM x JDIM x KDIM) ! dtemp Tendency of temperature due to dissipation of ke ! ! TAUSF = "CLIPPED" SAT MOMENTUM FLUX ( AT HALF LEVELS below top) -! +! !=================================================================== !----------------------------------------------------------------------- @@ -178,7 +190,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & !======================================================================= ! (Intent local) real , dimension(size(uwnd,1),size(uwnd,2)) :: xn, yn, psurf,ptop,taub - real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: theta + real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: theta real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)+1) :: taus real vsamp integer, dimension (size(uwnd,1),size(uwnd,2)) :: ktop, kbtm @@ -195,18 +207,18 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & ! -> -AETA(L)*PS*UMAG(L)*D(L)*GMAX/XL ! THETA POTENTIAL temperature at full model levels ! (dimensioned IDIM x JDIM x KDIM) -! PSURF Surface pressure +! PSURF Surface pressure ! (dimensioned IDIM x JDIM) ! PTOP Pressure at top of low-level layer ! (dimensioned IDIM x JDIM) ! KTOP Top model level index included in low-level layer ! (dimensioned IDIM x JDIM) ! KBTM Bottom model level index included in low-level layer -! usually the lowest level +! usually the lowest level ! (dimensioned IDIM x JDIM) !----------------------------------------------------------------------- ! type loop indicies - integer i, j, k, kd, kb, kt, kbp1, ktm1 + integer i, j, k, kd, kb, kt, kbp1, ktm1 !----------------------------------------------------------------------- ! Local variables needed only for code that ! implements supersource-like gravity wave drag. @@ -217,7 +229,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & real, dimension(size(uwnd,1),size(uwnd,2)) :: ulow, vlow, tlow, thlow real, dimension(size(uwnd,1),size(uwnd,2)) :: rlow, zsvar, bvfreq, x real, dimension(size(uwnd,1),size(uwnd,2)) :: depth, ave_p -integer, dimension(size(uwnd,1),size(uwnd,2)) :: ntop +integer, dimension(size(uwnd,1),size(uwnd,2)) :: ntop real, dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: th, sh_ang, test real, dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: sigma, del_sigma !real, dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)+1) :: sigma_half @@ -229,7 +241,6 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & kdim = size( uwnd, 3 ) kdimm1 = kdim - 1 kdimp1 = kdim + 1 - !----------------------------------------------------------------------- ! CODE VARIABLES DESCRIPTION @@ -259,7 +270,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & ! BNV,BNVK = "LOW-LEVEL",V. PROFILE - BRUNT VAISALA FREQ(1 ! (= N,N ! BNV2,BNVK2 = N**2, N(L)**2 -! HPRIME = Sub-grid scale mountain height +! HPRIME = Sub-grid scale mountain height ! over local domain (IDIM x JDIM) ! XL = EFFECTIVE MOUNTAIN LENGTH = (100KM EVERYWHERE) ! SIGTOP = HIGHEST LEVEL TO WHICH GRAVITY WAVE @@ -305,8 +316,8 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & ! (input via namelist), find highest model level. ptop(:,:) = (1.-low_lev_frac)*psurf(:,:) - do kd=kdim,1,-1 - where (pfull(:,:,kd) .ge. ptop(:,:)) + do kd=kdim,1,-1 + where (pfull(:,:,kd) .ge. ptop(:,:)) ktop(:,:) = kd end where end do @@ -328,7 +339,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & & xn,yn,taub,pfull, phalf,zfull,zhalf,vsamp,taus) ! calculate mountain gravity wave drag tendency contributions - call mgwd_tend (is,js,xn,yn,taub,phalf,taus,dtaux,dtauy, tausf) + call mgwd_tend (is,js,xn,yn,taub,pfull,phalf,taus,dtaux,dtauy,tausf) else if ( do_mcm_mg_drag ) then @@ -378,9 +389,9 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & do j = 1,jdim do i = 1,idim do k = ntop(i,j), kdim - ulow(i,j) = ulow(i,j) + del_sigma(i,j,k)* uwnd(i,j,k) - vlow(i,j) = vlow(i,j) + del_sigma(i,j,k)* vwnd(i,j,k) - tlow(i,j) = tlow(i,j) + del_sigma(i,j,k)* temp(i,j,k) + ulow(i,j) = ulow(i,j) + del_sigma(i,j,k)* uwnd(i,j,k) + vlow(i,j) = vlow(i,j) + del_sigma(i,j,k)* vwnd(i,j,k) + tlow(i,j) = tlow(i,j) + del_sigma(i,j,k)* temp(i,j,k) thlow(i,j) = thlow(i,j) + del_sigma(i,j,k)*th(i,j,k) depth(i,j) = depth(i,j) + del_sigma(i,j,k) end do @@ -440,7 +451,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & 20 continue if ( klast /= 0 ) test(i,j,1:klast) = 0.0 kcrit = klast + 1 - x(i,j) = x(i,j)/(1.-sigma(i,j,kcrit)) + x(i,j) = x(i,j)/(1.-sigma(i,j,kcrit)) ! should be x(i,j) = x(i,j)/(1 - sigma_half(klast)/sigma_half(kdim)) end do end do @@ -464,15 +475,20 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & dtemp = 0.0 endif +! Temperature tendencies crashing model due to size? Add a minimum size at x10-11. +where (dtemp<=1E-11) + dtemp = 0.0 +end where + return end subroutine mg_drag !======================================================================= -!############################################################################# - +!############################################################################# + subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & theta,xn,yn,taub) - + !------------------------------------------------------------------- @@ -493,12 +509,12 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & !======================================================================= ! (Intent local) real , dimension(size(uwnd,1),size(uwnd,2)) :: sumw, delp, ulow, bnv, & - & hprime, fr, g, ubar, vbar, bnv2 + & hprime, fr, g, ubar, vbar, bnv2 real grav2, xli, a, small integer idim, jdim,kdim,ie, je !----------------------------------------------------------------------- ! type loop indicies - integer i, j, k, kb, kt, kbp1, ktm1 + integer i, j, k, kb, kt, kbp1, ktm1 !----------------------------------------------------------------------- !=================================================================== @@ -516,7 +532,7 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & ! define local scalar variables xli=1.0/xl_mtn grav2=grav*grav - a = acoef + a = acoef !----------------------------------------------------------------------- @@ -565,7 +581,7 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & ! v197 uses p* as reference vlues for theta, in above 1000 hPa is used ! theta(:,:,:)=temp(:,:,:)*(pfull(:,:,:)/ & ! & phalf(:,:,kdim+1))**(-kappa) - + do j=1,jdim do i=1,idim kt=ktop(i,j) @@ -582,7 +598,7 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & ! ----------------------------------------- small = epsilon(ulow) - where (bnv2(:,:) .gt. 0.0) + where (bnv2(:,:) .gt. 0.0) bnv(:,:) = sqrt(bnv2(:,:)) fr (:,:) = bnv(:,:)*hprime(:,:)/(ulow(:,:) + small) g (:,:) = gmax*fr(:,:)*fr(:,:)/(fr(:,:)*fr(:,:)+a*a) @@ -596,7 +612,7 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & end subroutine mgwd_base_flux -!############################################################################# +!############################################################################# subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & xn,yn,taub,pfull,phalf,zfull,zhalf,vsamp,taus) @@ -606,7 +622,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & real, intent(in), dimension (:,:,:) :: & & uwnd, vwnd, temp, theta, pfull, phalf,zfull, zhalf real, intent(in), dimension (:,:) :: xn, yn, taub - real, intent(in) :: vsamp + real, intent(in) :: vsamp integer, intent(in), dimension (:,:) :: ktop, kbtm !=================================================================== ! Arguments (intent out) @@ -614,16 +630,16 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & !======================================================================= ! (Intent local) real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: & - & dterm, dudz + & dterm, dudz real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)+1) :: & & umag, bnvk2, d,d2, d2i, d2udz2, extend real grav2, xli, small integer :: idim, jdim, kdim, kdimm1, kdimp1 !----------------------------------------------------------------------- ! type loop indicies - integer i, j, k, kb, kt, kbp1, ktm1 + integer i, j, k, kb, kt, kbp1, ktm1 !----------------------------------------------------------------------- -! type flux cutoff +! type flux cutoff integer kcut !======================================================================= @@ -740,7 +756,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! compute wkb extension term for umag > 0 ! --------------------------------------- - where (umag(:,:,:).gt.0.0) + where (umag(:,:,:).gt.0.0) extend(:,:,:) = vsamp*d2udz2(:,:,:)/umag(:,:,:) elsewhere extend(:,:,:) = 0.0 @@ -784,7 +800,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! initialize d2i to a large number, which will result in a very ! small vertical wavelength (d) where umag = 0. - where (umag(:,:,:).gt.0.0) + where (umag(:,:,:).gt.0.0) d2i(:,:,:) = (bnvk2(:,:,:)/(umag(:,:,:)* & & umag(:,:,:)) - extend(:,:,:) ) elsewhere @@ -797,7 +813,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! for 1/d**2 approaching 0 calculate d by dividing by a ! very small but finite number - where (d2i(:,:,:) .lt. 1.e-30) + where (d2i(:,:,:) .lt. 1.e-30) d(:,:,:) = 1.e+30 elsewhere d2(:,:,:) = 1./d2i(:,:,:) @@ -806,7 +822,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! set d=0 for umag=0. - where (umag(:,:,:).eq.0.0) + where (umag(:,:,:).eq.0.0) d(:,:,:) = 0.0 endwhere @@ -840,16 +856,16 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! keep taus profile constant across top model layer (original code) ! calculate taus profile in top model layer (new code) - + taus(:,:,1) = taus(:,:,2) !del taus(:,:,1) = -phalf(:,:,1)*umag(:,:,1)*umag(:,:,1) & !del & *d(:,:,1)*xli*gmax / (temp(:,:,1)*rdgas) ! do not allow wave breaking for unstable layers - + do k = 1,kdimp1 - where ( bnvk2(:,:,k) .lt. 0.0) + where ( bnvk2(:,:,k) .lt. 0.0) taus(:,:,k) = taub(:,:) endwhere end do @@ -858,8 +874,8 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! ------------------------------------------- ! tausat(:,:,1) = 0. ! use all forcing ! Instead, let remaining flux escape above flux_cut_level - if( flux_cut_level > 0.0 ) then - kcut= 1 + if( flux_cut_level > 0.0 ) then + kcut= 1 do while( phalf(1,1,kcut) < flux_cut_level ) kcut= kcut+1 enddo @@ -871,13 +887,13 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & end subroutine mgwd_satur_flux -!############################################################################# +!############################################################################# -subroutine mgwd_tend (is,js,xn,yn,taub,phalf,taus,dtaux,dtauy,tausf) +subroutine mgwd_tend (is,js,xn,yn,taub,pfull,phalf,taus,dtaux,dtauy,tausf) !=================================================================== ! Arguments (intent in) - real, intent(in), dimension (:,:,:) :: phalf, taus + real, intent(in), dimension (:,:,:) :: pfull,phalf, taus real, intent(in), dimension (:,:) :: xn, yn, taub integer, intent(in) :: is, js !=================================================================== @@ -885,9 +901,15 @@ subroutine mgwd_tend (is,js,xn,yn,taub,phalf,taus,dtaux,dtauy,tausf) real, intent(out), dimension (:,:,:) :: dtaux, dtauy, tausf !======================================================================= ! (Intent local) - real , dimension(size(phalf,1),size(phalf,2),size(phalf,3)) :: dterm + real , dimension(size(phalf,1),size(phalf,2),size(phalf,3)) :: dterm real , dimension(size(phalf,1),size(phalf,2),size(phalf,3)+1) :: taup integer kdim, kdimp1 + + ! (Intent local) (trop/strat calc) + real , dimension(size(dtaux,1),size(dtaux,2)) :: psurf,ptop_tp + integer, dimension (size(dtaux,1),size(dtaux,2)) :: kbtm + real pfull_av, ptop_av + integer ktop_tp !----------------------------------------------------------------------- ! type loop indicies integer k, kd @@ -919,32 +941,128 @@ subroutine mgwd_tend (is,js,xn,yn,taub,phalf,taus,dtaux,dtauy,tausf) ! <><><><><><><><><><><><><><><><><><><><><><><><><><><><> ! <><><><><><><><> DE-CELERATION CODE <><><><><><><><> - + ! CALCULATE DECELERATION TERMS - DTAUX,DTAUY ! ------------------------------------------ do k=1,kdim dterm(:,:,k) = grav*(taup (:,:,k+1)-taup (:,:,k)) & & /(phalf(:,:,k+1)-phalf(:,:,k)) - + dtaux(:,:,k) = xn(:,:)*dterm(:,:,k) dtauy(:,:,k) = yn(:,:)*dterm(:,:,k) end do - + ! ! RC Added to make stable at higher values of gmax (above 1) REMOVE AS UNREALISTIC!! + ! where (dtaux>0.001) + ! dtaux = 0.001 + ! end where + + ! where (dtaux<-0.001) + ! dtaux = -0.001 + ! end where ! print sample output ! print*, ' mgdrag output for i,j=', is,js -! print *,'taub = ', taub(is,js) -! print *,'taus = ', taus(is,js,:) -! print *,'taup = ', taup(is,js,:) +! print *,'taub = ', taub(is,js) +! print *,'taus = ', taus(is,js,:) +! print *,'taup = ', taup(is,js,:) -! *********************************************************** +!----------------------------------------------------------------------- +!! Turning off drag for stratosphere/troposphere +!----------------------------------------------------------------------- +! calculate bottom of low-leveles layer = lowest level unless kbot is present +! copied this code from the low level fraction calculation in mg_drag + ! if (present(kbot)) then + ! kbtm(:,:) = kbot(:,:) + ! else + ! kbtm(:,:) = kdim + ! endif + ! ! calculate top of low-level layer, first get surface p from phalf + ! if (present(kbot)) then + ! do j=1,jdim + ! do i=1,idim + ! psurf(i,j) = phalf(i,j,kbtm(i,j)+1) + ! end do + ! end do + ! else + ! psurf(:,:) = phalf(:,:,kdimp1) + ! endif + psurf(:,:)=pfull(:,:,kdim) +! Based on fraction of model atmosphere to be considered below the "tropopause" +! (input via namelist), find highest model level in the "troposphere". +! This definition is to do with the distribution of drag, not the actual tropopause + + ptop_tp(:,:) = (1.-tp_level)*psurf(:,:) + ptop_av=sum(ptop_tp(:,:))/size(ptop_tp(:,:)) + do kd=kdim,1,-1 + pfull_av=sum(pfull(:,:,kd))/size(pfull(:,:,kd)) + if (pfull_av .ge. ptop_av) then + ktop_tp = kd + endif + end do +! Turn off drag below the "tropopause" + if (no_trop_drag) then + do k=kdim,ktop_tp,-1 + dtaux(:,:,k) = 0 + dtauy(:,:,k) = 0 + end do + endif + +! Turn off drag above the "tropopause" + if (no_strat_drag) then + do k=ktop_tp,1,-1 + dtaux(:,:,k) = 0 + dtauy(:,:,k) = 0 + end do + endif + +!! Idealised Drag + if (do_idealised_NH) then + dtaux(:,:,:) = 0 + dtauy(:,:,:) = 0 + + if (ideal_lat_cent_NH .GE. js .AND. ideal_lat_cent_NH .LE. js+size(dtaux, 2 )-1) then + dtaux(:, MOD(ideal_lat_cent_NH-1,size(dtaux,2))+1, ideal_pfull_cent) = ideal_strength + dtaux(:, MOD(ideal_lat_cent_NH-1,size(dtaux,2))+1, ideal_pfull_cent-1) = ideal_strength + + if (MOD(ideal_lat_cent_NH,size(dtaux,2)) .NE. 0) then + dtaux(:, MOD(ideal_lat_cent_NH-1,size(dtaux,2))+2, ideal_pfull_cent) = ideal_strength + dtaux(:, MOD(ideal_lat_cent_NH-1,size(dtaux,2))+2, ideal_pfull_cent-1) = ideal_strength + endif + endif + + if (ideal_lat_cent_NH .EQ. js-1 ) then + dtaux(:,1,ideal_pfull_cent) = ideal_strength + dtaux(:,1,ideal_pfull_cent-1) = ideal_strength + endif + endif + + if (do_idealised_SH) then + dtaux(:,:,:) = 0 + dtauy(:,:,:) = 0 + + if (ideal_lat_cent_SH .GE. js .AND. ideal_lat_cent_SH .LE. js+size(dtaux, 2 )-1) then + dtaux(:, MOD(ideal_lat_cent_SH-1,size(dtaux,2))+1, ideal_pfull_cent) = ideal_strength + dtaux(:, MOD(ideal_lat_cent_SH-1,size(dtaux,2))+1, ideal_pfull_cent-1) = ideal_strength + + if (MOD(ideal_lat_cent_SH,size(dtaux,2)) .NE. 1) then + dtaux(:, MOD(ideal_lat_cent_SH-1,size(dtaux,2)), ideal_pfull_cent) = ideal_strength + dtaux(:, MOD(ideal_lat_cent_SH-1,size(dtaux,2)), ideal_pfull_cent-1) = ideal_strength + endif + endif + + if (ideal_lat_cent_SH .EQ. js+size(dtaux,2)) then + dtaux(:,size(dtaux,2),ideal_pfull_cent) = ideal_strength + dtaux(:,size(dtaux,2),ideal_pfull_cent-1) = ideal_strength + endif + endif +!----------------------------------------------------------------------- end subroutine mgwd_tend !####################################################################### - subroutine mg_drag_init( lonb, latb, hprime ) + subroutine mg_drag_init(is, ie, js, je, lonb, latb, hprime ) !======================================================================= ! ***** INITIALIZE Mountain Gravity Wave Drag @@ -955,14 +1073,15 @@ subroutine mg_drag_init( lonb, latb, hprime ) ! lonb = longitude in radians of the grid box edges ! latb = latitude in radians of the grid box edges !--------------------------------------------------------------------- + integer, intent(in) :: is, ie, js, je real, intent(in), dimension(:) :: lonb, latb - + !--------------------------------------------------------------------- ! Arguments (Intent out - optional) ! hprime = array of sub-grid scale mountain heights !--------------------------------------------------------------------- real, intent(out), dimension(:,:), optional :: hprime - + !--------------------------------------------------------------------- ! (Intent local) !--------------------------------------------------------------------- @@ -983,7 +1102,7 @@ subroutine mg_drag_init( lonb, latb, hprime ) unit = open_namelist_file() ierr = 1 do while( ierr .ne. 0 ) - read ( unit, nml = mg_drag_nml, iostat = io, end = 10 ) + read ( unit, nml = mg_drag_nml, iostat = io, end = 10 ) ierr = check_nml_error(io,'mg_drag_nml') end do 10 continue @@ -1007,8 +1126,8 @@ subroutine mg_drag_init( lonb, latb, hprime ) ix = size(lonb(:)) - 1 iy = size(latb(:)) - 1 - allocate( Ghprime(ix,iy) ) ; Ghprime = 0.0 - + allocate( Ghprime(is:ie,js:je) ) ; Ghprime = 0.0 + !------------------------------------------------------------------- module_is_initialized = .true. !--------------------------------------------------------------------- @@ -1041,9 +1160,13 @@ subroutine mg_drag_init( lonb, latb, hprime ) ' is not a valid value for source_of_sgsmtn', FATAL) endif + ! ! RC Added to make stable at higher values of gmax (above 1) 1500 STILL HAD ISSUES, TRY 1200 + ! where (Ghprime>700) + ! Ghprime = 700 + ! end where ! return sub-grid scale topography? if (present(hprime)) hprime = Ghprime - + !===================================================================== end subroutine mg_drag_init diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 81963dbc6..015f180a3 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -627,7 +627,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l if(do_damping) then call pressure_variables(p_half_1d,ln_p_half_1d,pref(1:num_levels),ln_p_full_1d,PSTD_MKS) pref(num_levels+1) = PSTD_MKS - call damping_driver_init(rad_lonb_2d(:,1),rad_latb_2d(1,:), pref(:), get_axis_id(), Time, & + call damping_driver_init(is, ie, js, je, rad_lonb_2d(:,1),rad_latb_2d(1,:), pref(:), get_axis_id(), Time, & sgsmtn) !note that in the original this is pref(:,1), which is the full model pressure levels and !the surface pressure at the bottom. There is pref(:2) in this version with 81060 as From a1d4b05ff392821851e114e3ccde5fc3e7499177 Mon Sep 17 00:00:00 2001 From: sit23 Date: Thu, 13 Aug 2026 14:50:56 +0100 Subject: [PATCH 2/5] Fix mg_drag restart I/O writing/reading only one PE's local domain tile read_data/write_data for the Ghprime (subgrid-orography variance) restart field omitted the grid_domain argument that every other domain-decomposed Isca field passes (c.f. hs_forcing.F90's tg_prev). Without it, write_data at mg_drag_end doesn't gather all PEs' local latitude bands into the true global field -- it writes just one PE's local tile. The next month's read_data then loads that malformed file into every PE's local array (shapes coincidentally match under a uniform decomposition), so the global field becomes one PE's tile repeated across the latitude axis instead of the real topography. Only affects source_of_sgsmtn='input' runs chained across multiple restart months -- source_of_sgsmtn='computed' recomputes fresh each month rather than reading back the previous restart, so it never exercised this path. Verified with a 2-month smoke test: month 2 is now bit-identical to month 1 (previously diverged completely, correlation ~0). --- src/atmos_param/mg_drag/mg_drag.f90 | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/atmos_param/mg_drag/mg_drag.f90 b/src/atmos_param/mg_drag/mg_drag.f90 index 1e0ac7701..99f3fbef1 100644 --- a/src/atmos_param/mg_drag/mg_drag.f90 +++ b/src/atmos_param/mg_drag/mg_drag.f90 @@ -13,9 +13,11 @@ module mg_drag_mod use fms_mod, only: mpp_npes, field_size, file_exist, write_version_number, stdlog, & mpp_pe, mpp_root_pe, error_mesg, FATAL, NOTE, read_data, write_data, & - open_namelist_file, close_file, check_nml_error, open_restart_file, mpp_error + open_namelist_file, close_file, check_nml_error, open_restart_file, mpp_error, & + set_domain use fms_io_mod, only: get_restart_io_mode use constants_mod, only: Grav, Kappa, RDgas, cp_air + use transforms_mod, only: grid_domain !----------------------------------------------------------------------- implicit none @@ -1144,7 +1146,8 @@ subroutine mg_drag_init(is, ie, js, je, lonb, latb, hprime ) if ( file_exist( 'INPUT/mg_drag.res.nc' ) ) then if (mpp_pe() == mpp_root_pe()) call mpp_error ('mg_drag_mod', & 'Reading NetCDF formatted restart file: INPUT/mg_drag.res.nc', NOTE) - call read_data ('INPUT/mg_drag.res.nc', 'ghprime', Ghprime) + call set_domain(grid_domain) + call read_data ('INPUT/mg_drag.res.nc', 'ghprime', Ghprime, grid_domain) else if ( file_exist( 'INPUT/mg_drag.res' ) ) then if (mpp_pe() == mpp_root_pe()) call mpp_error ('mg_drag_mod', & 'Reading native formatted restart file.', NOTE) @@ -1179,7 +1182,8 @@ subroutine mg_drag_end if(do_netcdf_restart) then if (mpp_pe() == mpp_root_pe()) call mpp_error ('mg_drag_mod', & 'Writing NetCDF formatted restart file: RESTART/mg_drag.res.nc', NOTE) - call write_data('RESTART/mg_drag.res.nc', 'ghprime', ghprime) + call set_domain(grid_domain) + call write_data('RESTART/mg_drag.res.nc', 'ghprime', ghprime, grid_domain) else if (mpp_pe() == mpp_root_pe()) call mpp_error ('mg_drag_mod', & 'Writing native formatted restart file.', NOTE) From ddc876832cff83a009b410073f55d2c90834dd9e Mon Sep 17 00:00:00 2001 From: sit23 Date: Thu, 13 Aug 2026 15:24:32 +0100 Subject: [PATCH 3/5] Add cg_drag/mg_drag test cases, trip test registration, and ERA5 regrid script Moves the two GWD on/off comparison test cases (previously ad-hoc scripts in isca-experiments) into exp/test_cases/ proper, restructured to match the repo's standard single-exp convention (module-level exp with the scheme on, cb.compile() deferred to __main__) so they work with the trip test harness: - exp/test_cases/cg_drag/cg_drag_qbo_test_case.py - exp/test_cases/mg_drag/mg_drag_socrates_topo_test_case.py + its input/mg_drag.res.nc (ERA5-derived subgrid-orography stdev, T42, ascending latitude -- see the regrid script's docstring) Registers both with exp/test_cases/trip_test/trip_test_functions.py (get_nml_diag + list_all_test_cases_implemented_in_trip_test), verified by importing both test cases through get_nml_diag directly. Also relocates the ERA5 sdor -> T42 regridding script from isca-experiments to src/extra/python/scripts/regrid_era5_sdor_to_t42.py, updating its output path to write directly into the new mg_drag test case's input/ directory. --- .../cg_drag/cg_drag_qbo_test_case.py | 192 +++++++++++++++ exp/test_cases/mg_drag/input/mg_drag.res.nc | Bin 0 -> 74950 bytes .../mg_drag_socrates_topo_test_case.py | 228 ++++++++++++++++++ .../trip_test/trip_test_functions.py | 16 ++ .../python/scripts/regrid_era5_sdor_to_t42.py | 111 +++++++++ 5 files changed, 547 insertions(+) create mode 100644 exp/test_cases/cg_drag/cg_drag_qbo_test_case.py create mode 100644 exp/test_cases/mg_drag/input/mg_drag.res.nc create mode 100644 exp/test_cases/mg_drag/mg_drag_socrates_topo_test_case.py create mode 100644 src/extra/python/scripts/regrid_era5_sdor_to_t42.py diff --git a/exp/test_cases/cg_drag/cg_drag_qbo_test_case.py b/exp/test_cases/cg_drag/cg_drag_qbo_test_case.py new file mode 100644 index 000000000..fc9a93403 --- /dev/null +++ b/exp/test_cases/cg_drag/cg_drag_qbo_test_case.py @@ -0,0 +1,192 @@ +""" +Non-orographic (convective) gravity wave drag test case, using the cg_drag +scheme -- originally implemented by Chaim Garfinkel (based on Alexander & +Dunkerton, JAS 1999), with edits by Stephen Thomson and further development +by Ross Castle (see src/atmos_param/cg_drag/cg_drag.f90 for the parameter +history). + +Physics settings (do_rayleigh off, do_cg_drag on, tuned cg_drag_nml) are the +QBO-capable configuration: do_rayleigh is off so the two forcings don't fight +each other. NUM_MONTHS is a short comparison run -- long enough for the +cg_drag forcing to develop and show up in the zonal-mean wind, well short of +the multi-year integration a full QBO cycle would need. +""" +import os + +from isca import IscaCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE + +NCORES = 16 +RESOLUTION = 'T42', 50 +NUM_MONTHS = 6 + +cb = IscaCodeBase.from_directory(GFDL_BASE) + +exp = Experiment('cg_drag_qbo_test', codebase=cb) +exp.clear_rundir() + +exp.inputfiles = [os.path.join(GFDL_BASE, 'input/rrtm_input_files/ozone_1990.nc')] + +diag = DiagTable() +diag.add_file('atmos_monthly', 30, 'days', time_units='days') +diag.add_field('dynamics', 'ps', time_avg=True) +diag.add_field('dynamics', 'bk') +diag.add_field('dynamics', 'pk') +diag.add_field('atmosphere', 'precipitation', time_avg=True) +diag.add_field('mixed_layer', 't_surf', time_avg=True) +diag.add_field('dynamics', 'sphum', time_avg=True) +diag.add_field('dynamics', 'ucomp', time_avg=True) +diag.add_field('dynamics', 'vcomp', time_avg=True) +diag.add_field('dynamics', 'temp', time_avg=True) +diag.add_field('rrtm_radiation', 'co2', time_avg=True) +diag.add_field('damping', 'udt_cgwd', time_avg=True) # cg_drag zonal wind tendency +exp.diag_table = diag + +exp.namelist = namelist = Namelist({ + 'main_nml': { + 'days': 30, + 'hours': 0, + 'minutes': 0, + 'seconds': 0, + 'dt_atmos': 120, + 'current_date': [1, 1, 1, 0, 0, 0], + 'calendar': 'thirty_day' + }, + + 'idealized_moist_phys_nml': { + 'do_damping': True, + 'turb': True, + 'mixed_layer_bc': True, + 'do_virtual': False, + 'do_simple': True, + 'roughness_mom': 3.21e-05, + 'roughness_heat': 3.21e-05, + 'roughness_moist': 3.21e-05, + 'two_stream_gray': False, # Use RRTM, not grey radiation + 'do_rrtm_radiation': True, + 'convection_scheme': 'FULL_BETTS_MILLER' + }, + + 'vert_turb_driver_nml': { + 'do_mellor_yamada': False, + 'do_diffusivity': True, + 'do_simple': True, + 'constant_gust': 0.0, + 'use_tau': False + }, + + 'diffusivity_nml': { + 'do_entrain': False, + 'do_simple': True, + }, + + 'surface_flux_nml': { + 'use_virtual_temp': False, + 'do_simple': True, + 'old_dtaudv': True + }, + + 'atmosphere_nml': { + 'idealized_moist_model': True + }, + + 'mixed_layer_nml': { + 'tconst': 285., + 'prescribe_initial_dist': True, + 'evaporation': True, + 'depth': 100., + 'albedo_value': 0.25, + 'do_qflux': False, + }, + + 'betts_miller_nml': { + 'rhbm': .7, + 'do_simp': False, + 'do_shallower': True, + }, + + 'lscale_cond_nml': { + 'do_simple': True, + 'do_evap': True + }, + + 'sat_vapor_pres_nml': { + 'do_simple': True + }, + + 'damping_driver_nml': { + 'do_rayleigh': False, # off, so cg_drag is the only source of GWD forcing + 'trayfric': -0.5, + 'sponge_pbottom': 50., + 'do_conserve_energy': True, + 'do_mg_drag': False, + 'do_cg_drag': True, + }, + + 'cg_drag_nml': { + 'Bt_0': 0.0043, + 'Bt_nh': 0.0014, + 'Bt_eq': 0.0043, + 'Bt_sh': 0.0014, + 'phi0n': 15., + 'phi0s': -15., + 'dphin': 10., + 'dphis': -10., + 'flag': 1, + 'Bw': 0.4, + 'Bn': 8.4, + 'cw': 35.0, + 'cwtropics': 35.0, + 'cn': 2.0, + 'kelvin_kludge': 1.0, + 'weighttop': 0.7, + 'weightminus1': 0.28, + 'weightminus2': 0.02, + 'source_level_pressure': 315.e+02, + 'damp_level_pressure': 0.85e+02, + 'cg_drag_freq': 21600 + }, + + 'rrtm_radiation_nml': { + 'do_read_ozone': True, + 'ozone_file': 'ozone_1990', + 'solr_cnst': 1360., + 'dt_rad': 3600, + 'do_read_co2': False, + 'co2ppmv': 300, + }, + + 'diag_manager_nml': { + 'mix_snapshot_average_fields': False + }, + + 'fms_nml': { + 'domains_stack_size': 600000 + }, + + 'fms_io_nml': { + 'threading_write': 'single', + 'fileset_write': 'single', + }, + + 'spectral_dynamics_nml': { + 'damping_order': 4, + 'water_correction_limit': 200.e2, + 'reference_sea_level_press': 1.0e5, + 'num_levels': 50, + 'valid_range_t': [100., 800.], + 'initial_sphum': [2.e-6], + 'vert_coord_option': 'uneven_sigma', + 'surf_res': 0.1, + 'scale_heights': 7.9, + 'exponent': 1.4, + 'robert_coeff': 0.03, + } +}) + +exp.set_resolution(*RESOLUTION) + +if __name__ == '__main__': + cb.compile() + exp.run(1, use_restart=False, num_cores=NCORES) + for i in range(2, NUM_MONTHS + 1): + exp.run(i, num_cores=NCORES) diff --git a/exp/test_cases/mg_drag/input/mg_drag.res.nc b/exp/test_cases/mg_drag/input/mg_drag.res.nc new file mode 100644 index 0000000000000000000000000000000000000000..dbaf842f77d5a114b7942738ab174b05431e7539 GIT binary patch literal 74950 zcmeEv2Ur!!()J-EIp>^n&S{FAbIw_E4uSy?6af)Y5d~3^V8%ca5lIpZfT)-N5s)O3 zM6w9~VL9HtS3jTq_P={~-S0Bb1MjKs>guYudgk=ZndwQ@)6yX&Vj{vOBt(7P)-hWP z57sY55}1lK+l&ZPJuM4tl5%XWU&?TBF;!oie2YZbs9$VQP;4kZCEZll!nE21S{BrX z85}Isx2?a*o2_5rZ0cKLTI17ge89rR#>hB`m8QM2=`B&ex{DSA0t0Q_1EakE zM;bRs{~E(r+_bY4k0u$nv3xU0Z2y#`giQt1SN#7;k}%1pky993e@c=z0{r8X6cUU! zs}n!5@B4ggbSO@7ro&WI&yYlWBh-yV{xUjTO*GBILqfg7BmL0BYh#pfq5n3r9_<|- z;TICDDkT#9%Q#e(5tEbDI3yJ`t(;+VvKbolotQwdy@_NZR>dXw-}HUI||aRj-%I5_{4_4E4qxJ4Shp zRqG#S8MZbgD%BS*6q-9%$=3&|kYz1R%G6g1TjU<&ma0D{QSip8RIEOitw{DLt6+UR zVS;gM3U~d(v(bJ_qs;aEk6S*u*HY9^7rwC(wOpvfAc83ej2>Xx0mcqs>;k4A!1N0k zKLFzwVEhD(-++k&FmVAUPQb(s7#;w_3t)Hx3~zwp5iq<0hG)R=4wyUuCNF@=6JYWN zm^=a|uYk!jVDb*=f%5B?fDgY233Qv-&w-KoHjzUFQ!cTIH>5H3a+_$TxQQGpn>eGk ziT$vNdK#N}Lvs@uv^Oz9XA{SCH_=Rg6RQk1k;7;cQ;aup#$*#k%r>zfgIVUAvYy2z zu43?p<)-XqwTTSYo7jfI1e;A+#OFqmhzDVy1EBCf+GR$(yEaZ~1S+QcplrZ{iP zDlVHigTW%#OgScwr`dvV6Yp5GCrH-85kVHpoZ^ec^(EAFlgqtS$-LVxc-~63kIt&NENVI z9*Dsf3~~f+mhZ%17Y4I2Fozm zk3rI~O+6(H`d}~%gAEv*!XR__roJ8qqcB*2!N(X}#h_rsroIgZlQ4J#gKscM6uGG< zk3laC9>HKO20vhsA!<`!8-w8(Jd43L3@&4kFM3no5`(C%?Z1_QjVZ~|LYke7P`?NIb&=+|MBGJj1h%}#k1iNHhPK)8z0bH!Y}_(D8S%X`S;Br38^-k z<6yP?Vt?cFMieO2*EuBKk0<)t8+yNMqL01N0=<33jn5n1ZE{aXwjrZXa+@BYPz0Ou zFaNdM+2%IV^!3!}YZVT9UO$Zy@R>F~VBz9mWC9=G&~U#X?-r8H0^HB6@nvZL+dyw}pb5Lb;uN#p@M0xmx`+12*c)ACAbA*J4_=LNM z`o?lZM7jrixrcjkczH+rxksX>sT=~^!$X2Nv@O(SIU>A5!UaiKHisYPEGR12FLHB! zyFt;0?I04T%`+5o^obG8Kc5vD>04T%7||*e=s$Y&A0hg0;|%31`CoCvR13Q4#Pw?T72j*{5$z(QL`+d|t!Q`_|y59H8V8ym!c8?) zJ!WEq5w~EAfh`8M7}#Q9i-9c$wiwuAV2gn*2DTX3VqlAbEe3uf2C|HUgGpsSBcZP6 z`fl^hA*@3=y_a|Q4jg!%1YZ{nSwFqKp(~lLN!JVM>@A z5p*8NHld;-g_EUg?z@XgK#ch^{!T1BkO^bhT`9SQR4z|4m@>^EF`pet5OlXBZi(TXelNHHOeoortXL`-IHs_&rn4 zSx55xEVgxf;sKw*2&JnxA+Q~=Mb&f?z)4C;Hey~JKs6r7-m!lLX^{~>n-h9aoB#H@YCaqC#>Dl0sRe$%Rcs7usdl0*G7Yn4~IOyfo(gSI$&%HgKG+;`_ z?oKLC3uB*md^ktR0CSc?<>IH=U_x*r{!Wx2kJw3G*n-5g97e+^loBixbxx*+deZ&SaT(f6xg{9 z%DmcEZ8C5n^R@9sY6To%#Cz$Mdmax2uMA8+swM_{l27FX4bxpC(cVobBl-(R(W)ZXQ?^ndsTVGI8()nSec-m=$k;xXfb zzh9F9OAZQj4@U~0?*4+vU0aoCU;c!QUj9fP`{g52awNlcjNl`(q+>Aksc#P9$tgde zQ#yf|R6WIii#vrBYN6x>P7WesdF8@>8Y9S)Z6iVK2JaE6y+_-OpHCo)q7>L`I5UWM z4&TzNiVw(6)sUvg=Oz%jy6Y>dpQaGPakd1p*E5K3Z@{h3wkwFq{$r!1FBTCQZ8>I9 z+fRs{wXe$kwvWi`^;NYilq-k<2bou>EFJ{-pqjgG-~(Q40?(D3tH_B!$_VKbACa9j zu5dhi2EilVKTvgb4!Qm!A}s2ntpT6EAmcA=N$!t*Mk?+;b-7nFfQTfj4#p0QA{JV4 zPXZTvkuud9Q$5;#geKq+F7DfRNXwnbm22(&h|A^EuR5K3kfA1%frpMGNT_54?J7DB zX+JlXz<^K4^_PR?Va?MBA%T3^t|K21U+jeWoD0*)`h!C1rs63?eylFe?#(#DBb;F^ zd3+iPJ-Ucuy1#Q?$~A*SL4jg%0b$xN>3j zg?Lry&OAND#t?=tCPydP=!05wefqN`6L2^zT)d2J4D&+o z9D0qlLE-Iv{!`agpvT9vIeQI#Ym}(b)pI3EkcKnvDaoV&$Jt8UeCD*_scUhjUbPu; zr6as}%7*Y<{EIH}Co8BI+M+XHs{Hq^GyTEYG-f)lqMn1XSl#=U*>N^p&$KkdUM zZXi%h$jfl02a@_4YVImV=$<^$l5k-4O-8uzQQM)j*$E3YLf!3>y62I!kvFr?73L9#9ohly+zUv8 zBuCWcgc&5O&ZGDFfhj~n+NYUFb_U72c4@K_Zw7hKUB*z{HIF!HFSf6oT}1pT)+et# znL$QfNvu}Fut7rOxw}_4AFE4 zh4ge6KQ1&&P>5kpTs7VtWKtOnRthwI`($u@wD05^Sxln#Z*1CEckst7W)tott!0b&s!(AT`sl`U!YRfs~ z4#O?o?$0xbNoa9~CiXn?oO74>^UN8fGf@;&6^TIY%%kMg=eXdI#_rW{;Uj{lQ2fC6 z8Ww0B66k2|qX$<(zVA9x67it1=-3xz;p&^C8r(|=qbyQK1u-zd8qiSg9eZ!f2LT>))sQMtEydj5(b>Dr7siuc0j(og|%3w z0~km>Y7=_6ii9O+pG2S^q39(edH-Y((RA^PBFlM)u$|EpJTEbTEL@UhW=|SL8uQ;P z7v-!V17>Wm^i;X$jkqg51_73#%3c(%NiKoFJ z1ugV5G%<86Fx1~Pv)ID~EqV|6t*My6{dUkc{X|aa8o|zVDCP#c+V@X!edwT--K^Y7 zh74k{&XQn9k->^S`-5v$T=3*_`?-KPJ6M|qcUJwhSR%VF9tkc!b=1k*da9IR!5XG2UuS%$T=Xw0dAJ$Oa)z> zFiYoG)x${zRk3bwnpGDOjk_)elYK-GtvDxMM7E4*WS-qYX|#-lKl>6RK)r&rm*0JR z(v<*)NzJ8~?{UB*x#~C>EOt0jo9FXco*nM3E@rlnQ9?lfYc_7#9S{Fp;EIN@wniam)J zGwg0kIEW35(1_2bSkcN4UcB$z^Sju=U;4_kls-my;@}g1+Y=u;=#T?@IkmlU2NrUv^*SCHnwhZ)>n zL&)7BMVZn5=SXXY*dF@6Cy2XvED6h*SBPl+jEW%J8$<=?BBhspKN6aPrz^<(4oPmu zqqH>XK{!lgoy+ZCAYEFe5qboZNQ+y11ihjrlrqR`(YFGWg~u_ewJsvtWb2Bb=HkL` zzDMVC`-CC+_({~q5iPK4l<14s)&cRu2Xkt6D}vg~FT?>2LcpVXuxY?f5NNvddj)n% zfZoE{cd1o8VCa5p|7ag85Rjo>@a*9Oai$E1tY8-K)rX0FY1E*dXYhDcl@e~Gls&&# zOAJf5gWY*AQ9)&W&5`JNIxvuqTvyEC1YCalE4(OnsLK{^A{k-_F0Zt5yb@Xfxij-g z`1FviLeIv4u*2Bes1vQb5Og0BKYd446C~CZIhka9f&cckXPJ#Qz;l!oSCGO4a$FAX z7(Af}lnKw+?w==s=&Zo_s#A&}>B(r@b%Fz?^jFP0HRcg*+b6M3f&4HL)whdTdp$MfU{K z5fMGDR=S29=$-m-s=DdGM?%t+)a4Fm}4@mY;C0;i}HQZY#nC+%0pDU6VLNFfS`g4rvJ0!HQo zw3X2Fe8QZ!_{&M#!131Zy%ox&uu#mi`_>>CD5#Tcw_~FQ3np9<$tfl{^?F}nR2UPi z_>sDkH1fbxr;}DIvz#E(va(vWBmfKC+*%iwxL`5JE}G(i5WH5&A(*P?2brtwbOb~^ z;Os&pP+=tu7Dp!g25aabZ)qv7g^dQ>dhBT3Rw?0?4l_O@J`=FyR6niwhzt0(```J> zEFj)-#piX0NJ0Ai4Y+xj24)-S?sYeiz-!4?3n^nV_((KQoqcs1Oo#|{tI*;=nsXMT z+5QD2zrbXcI%^%bt5>iDYmBX$@7E#|Kh7 zSMRiJD){0kG;`aS8b0ccYUv-N1a|SfGP(>3XxYAQ;CE~rXvm%}gVO{M!&K+-iHaCn z+J$!H&#=H?On7T>8$S>e$X#(O_lDhrZy!YZi@{%gi3va1lf|}E3-DKaBk$vTLUYng zv6P=aoVJ_d+g7LpGrm|^ozsfIUn|70{!s!B~l5T|t}^m`ct>8dW2zutfD?tWWPJo%=`|DGzO87%VfB&dVXwqbJKJyO7Jc2vQ} z4+nbj_by*reuupKOm`%1a1=?bzCWAnF^UK|K3Pwt8$$$4((;3RiQpYx){NO3qTijL zm#)!oTUGb~d${PC>n&Y)>oqdCQ%42H_P&x5;Z%m(7jD(uixLJCUuNa8SgZ=hUhojL8o81fLXytM&@%jN3%u z#7Ta~iz@0snfP{x;+Q6ApKn1mObf!)ws4A;Biz8mXVL#QMFJ)d5#x~QlfbP41@(4! z95??Nq>rS9KA~sS2cEMs;r8 zAf6K1v~P<$$TLAgEJAW(mS|bzk8hh)yIvm{f|GybLmgAfEXjZCoNaB;OX%v zXHTCofs0<^FUw2?fOp`!$fYUW-;A%3O#-~0~*Xa2BT?XWEP zIZZq7-6INbP*umCEzrQ_&TS;LE3|M-YLHbP%q^FY5o)mTvI^L5rfSzCS>l&Iyw(Deu$MKu0Q>QI0MH#0+;$;}_6_7oLQP$PO;J ze8ZIcaTf#NV3({?45H_!hs6at>Y2d5^eh7nCo9bJj7&Hx@x$JM+kc^*a{C>X zgeYxjnm+PWnNJ0n=|2ehyLdt#A&pQNVh#HDpSsVHcmqkiw1I1PX%kZK>Vf`! z?d?PAhEQw&bm9cB0*qf&d}-LF4EyMBunS+&fG?LGMp?ycKm+ONB%WS5_`n#qimfIM zNyX``q0XA{0*mZMM7qd#{m|#1pArnx1_?bHLN#_nAoHYTr8bj*b*{FfyDvz><%v{_ z(PB|3XRMvmbrS%CV~PpRiky&pee9$Z`ueWA&WQse*BQW){{xA}eO@?8xPajt{dLCIMrlH9q21mx{UJBim1_8akmUAe6xa6Rg(Oks>&KIg`bmn;g zhjr?`e77`Iawnkj;{`!Z#OI-CqY${#=sB^L%7B^@*_sx!D7f(lGfVA4AlSk_s!`h< zP)F>Wg~|<~Kw%^_kc1K7^eGmiS?TZM{ZHp}%l^MIFhaLWIMY=hTnZSRTnjycJ5%k( zqBZeWUjCJ)-z`1d*XEt+De$NBD8KIYB8{F7W}@<{WD_P4i?CFC{74fhuygw~=PCG| zKXD)-vS**rcjdo3;H?IKYXtt>@zs0~G!~n%E#LVn1 zJN!ZBw#;lXu*JX@16vGiF|ftJ76V%hY%#FKz+b=sQJZtq&SpRGYCIPe-2BI%??;xJ z7*mCX!(R~ER+n20Y%%Zy2JT+am2B$sgWTx`SDA8U_<khkG2>O0Qn5#$dtzh!2Nfq!8jYmQd7tp>B?{lC=w zK0&uOVLZ z^nK44?YA9Uy#7vJx70To(4Nb(roSr&R}?<(B19=eZczHd`tCqDWx(?QIiwC-JpHu{ zbnIscVNmn|PYUYxNiTo!J_Kq*#}y!cHo8_u6M?mbya@*(4N#yL45*h?{Xd&8>HWv- zw#R;dzx`jEoxi^Kbn)_3B?I~2o%g3zxy99j=ye`0Y4!B&rvh>=5w$Oi*f4&Y2H(nt z8V<%>fkal*-);9>M*3>vL=dX}t?#FAiTRzfkfEeTB7ZE<)a?(VB&LED#xF)e?2OZL%DSH~{Qnw{PKFpZ?z7nmiSFhOw87%AtLSe1_Ms5SG+o2~(!>JHC4=VAMwtPJ z;EPGnO<7PF%wInz^4rfV|Lg=h9~-Q8zqtChMtrVaJURWJnl7&91C->g<3&BXz zELV#v1Gx7TIJUkPgollUwsai4(8*l!VWmVAkl5(p7EMl29$z{r_K5=Oz1OOr60iaj zGwlcYJa*uW&k;X+84HdK3XK!_{QK7f|EReCTkV*R5j!p7#Y4cIA(MxcvC#R#F)b%n z7Hr;>YxXJ5BR9LFolg|8KreqoKV6joFxlW=3vtkg!kX)AsMDur`R|zr@D!{cus$OIY6))s z)Hve6@3Z@ybbtiV7vr#`m5RUzo~tD_vw~_QI34I7Z_4lf8nMG1r|2xP}6Z+AbG|%tcG5niuPzj`Yt#0dYOb_BhnW|LHtvl z_pMfQxLT2v+(Tscx6R*0ae9h27xX%-vExJ8e4Ox%!}`EDdVQ-U{L(i&={R9bdp+mc z6Eirzve-T0stGl?FSAMXS$-IA>IFqR$`r#N)_=?U&nnlX;}0oz3d0_x=E#HB8gQY1 zL`kkm8s;S;{M%&3A?DT94n2DXDF5_G57$x>jyjz*@FYjCTh-@tKBz+yR(59U^?G2# zwaEM!PAPWK(S84Ry&Ao4b{aN^3B4w~3$?LS$iTef`e!BKzsRpBGH01%=d8i9RppfA z;tm)Z8tFIC_6N1%v7|%kvq;8~34BS`NrXcYXR2OX4W?cXpI44UU|eGWae847aoS#= zgS7uX$MHqA{gRfM77$*Py=%4wfc2eZXiyO!To zKd&`tt(6D;(wnc`MI_-zL*rsiQ=&pNdim;dT<^YE%q!u4r2IE1f0?V&yE9G@;;w$$ zyILg=pOzA})k4K!I<{oiZA}8Mmp#_^6c+)pw)CJd0UlV|PL&wNsS4rm_%dJnQo-8> zjqXWTWpHR>Q73lif~fHix;yd(fbhhefBWvObq0QxAwXrMcVsU<6{rOc`9C&xfYo>J zGO0#(K!lcKIAaGsoXrRb=ucljXjUR7&)Mie+BJRYyGK>QHxG-de7`(;J#i5t!=^y{Q-;#o)l&v_^clr1>?_2W!%s@SE=iI$s52)Vm zDwpqU5Az|{-!$-sLXbr-9_Qc!axmIHAY0NNEOhIvtPE$6TMSbLR$f-H{n;V>+s_PO z*A<>!9u)e(VECvp*~{*r*HLzv^xT^y3fJjA;82+9f;3ZFpdF4msM6Qg2u#}n_CVQb?J-lB zT|sEaW|h7hci~o7C=fM62&kt0QgrbD-nv%pN1Sx;(x`&O)MKAR(mLRag)`m1O$e-B z8PRcge*Zfv`*B$IBi(v_lBotmsvHr4pRLjBWU4^J zc2!t(=rDdxDF3WK03QAUw72yKFU>+@@_oI<5}+<3j9NEOyXTO3*hoXbX7J#$iQc|M_(X z|7yN|tAC>d&T(meTtHayHr0y{;9LFK(KFX#W4XR7|5k5H{&N`UZWnXI)!PU9gI#GM z5BC5bH5Om(X<5)}WA=D!Nd*g|C(T@gd0_3#+~HX@O*q{&o0ll24lHs?9+6qn(8he; zYdV4-7;^?PLqbt-;Qi6tew;dxzCc&5Zm0oVhp+?ETu$q__E-9mbPMtvIuM^9lG)52s$2d-I5-0Yj}9Rpk}Z@6D&R~+s7+; z!Zy~XmkNuP@bKj8^lJ3`PyU>&{dK)iWR1zXrq3Xu24WP}#s7Z1V!ERbxYGCjkhh%~ z#FeWzeScX0*D>*bqi46S3u>%8cfplniwqtY8PMo4DeVXs0msj$v4}5;!$r%s!$oKK zfcDmo!V^Ri@ZeBcq7q6KQm$)?_7jRgT>aV7od%L%v7_vTYrG7ImYWGQqSviYjuM_- z?>B&33|ivH(ANjg1rn5)rD#C>t1Q2g7>%DMp9a0xLhXb#;C-eSv#kO8J#73gXXB;N zzh~C%O?miM0dOg>iSgytpmX1q(~Ri(PRiu7kJ{+-8lewU`&f&bT&Lg|`z%GbrtF(K6J`k2YIWDU07l4O#pKmHiNrIqo`qc(u z7BD{EJ<41|3%PtmcVeb#LH6;<>?j9W;IO;&{%p4lEF8Y{SnfWcZ~LNr2bayDYk}Tk zFTF18@?t3&nKOVY>eUn7MH-+ya@!(J(H^|B{l{;=F#gN)UTM8om+LMe5G&)29hJ6& zzTKf>_`Q0t5MbN7he{o+XGI$B)&r0!VN0eiiUOC%8QZ71j_^o>uwKeU3nJo59tY*h zLK|-be?EgH?75hWm9g;i&2N5doEHY|tspb`@%2d|t)A>4Z=xUOi?}GU%z^0p^LY~G z8_zNd%plF5R^WYeEQBlMXJ2IZgKEZSR7!E^`-r~TytH1urkC3OUHLa^TjC$YK$@hE z;f>BnnCP&jIe8uh+IwgsX50ASm5k2BQnC=Z7_u#XcH{x_a(WxDr;Oml^P+r8l^u30 z79Q|oVgnLYfn!S1GQiaob1+p}77*{wiuFP%cU z#JGfh&p_S4Is3OZz*$GdU&vw!G3R`49I=;#F+MZ9E>j=qsc(GH=kVkE0sf%0{o7^| zlY@+O`sqMsuO>d5KmgqQAXYN)PuI(Zy@omWQ$Ip_ST)|CP&WdaBh2`P$_D?)+xDyb zUJ<_!{fGRmL#n4_9^YU7^}dHpKUS9MY^)#Zf7^ab{&N^0EEm%s9Z!Os-FgB=WZ{(x&(!!L=`3;0-U7nJDr4QL4adx9u&?Qa-4=&__$@sDYj*&`Mg)nF0RVr57dh zQc$z^dgGEvFcgg04pY-ef>z{*nfL5}+`Rd#zkf7dd#QY@6MY=dXey$aprV7X_IW!s z8wpuW;LN){w;xQ~!_@Tol+y`5z{K+8(v*+|9F-I}QkLlP5B;z0w(#%80PD1CT)l-N zINZQ5K7Jt-9%fRl=lTL9>o7TFzF>#TryTf?Rx^T-4=Hi`x*G(nzLQ6NbcgBGgD-;5 zJ3$pYek?74D5MxUikwZ71h%lmv3NFId~AB>=5qjvd6g8X&fO+MG2o4AhWg@4YmP zAoe(I_Qi62sPU@5)7~ZvS#-IC4qc)^mk_(2I4=QT`~9m>?sn7lI%;-U^V>m-uP6>w zXRqS=IeY!4UpL>!N`7wk1&98C!PnBiey@LxBm2C^66cBP@VoT_KI5^YK(KMZT1U+UJ~mmlBMn7_~fzARJuL9{75-nQ+tP<1>TBU)|?y%yUU6SD0=YsRtUmj|srXSeK zQo$NcE=O$w3f|ZU90_6z2feGUEj;b4|uBZ$(f#2Qk+kR4N80SGE9aJ%)im+fSJtz$d%R-FTO z4e_k%o>GE0w}+_~Wwn55E~ehxMj8f6yKmD3D!_rqJZn#8B|%-^KFUN)1!@v4qZ@i{ zVV@QArM*iwz;*lh6Cy@^xV*!;Q+;|4Sv)izfhEZX8YV+D`1$Hk#S`;}_JtH=jWC^8 zemjZ8@%lbe?jeG^n)okiDI_88a(C=0i2&%6_~7GZ>jzear(P&~xB~8=C|B{a7wkIQ z=`V6w3D|D)x(Af}alBc-$+=!H$MC!NyDMJ)y!33}|4vcY`~|x!vhsiZ{5%oUTE+W_ z@&C2{pEH`bdifCpA1hZ=-+c7@;dt>B&Prdt6%CbM#D|Ylgn_d8n{}-;J1BoP^2Ce@ z{d=C%v7W2f$bQ)VNB#cMN!nWNIZ>5 z54G-gP)A?qK337R+e;Rj-tZZNO4{||+Z&|+LB0aF^uenEmK*%Al`+?J|XG1 z=kF(DTm9p#plsLP!|Nn@npQN01WG(zrfij4bmenLX{qgr(|5;(a^zhT3yZ8~2ntz+0 zL&60kx31U)=IBC8KO=csvkoi~Y24-PmISXh`F46PVR-Izj)mWs7tTE@)Uaw2fb7+K z?W^0>z-v^F;Os$u;FwEf8q*?z+nqFg%TzcJhaI+!Ga3&ZCmhqIuFoJ{?~giWqhz4( z?NVR=j44D{2_HIpK^&$juBp|XR)W^gFFa1!^S~ktch);6HCRY~IFCyI_4s-s8{6o9 zzbA}Gd;1E6*u%D^6QL|qzrDZy7N0Jcsjl1MDu8>}cqjLY1N@f$xAg*Eji%iA?jU*8 zL9C_G5H1DXIVp026>{ua)XzK%{7rl0yA5UAcN_ht-nad3$$t(5!q2YO44(Ig3Z*k` z4a0jN=Y9ywVs9i6RGC>C)_*_m`*Xbd_j{p_yueMlypBYp%w`TSbHcrzK0}gD1uzN} z-{*|*f`4d-#*`2jeAK?aV>sUsv=!c)JLRK~2Zur~w+8S+RIDN2HZ?9Fbyxf-I&TRH zx(;}K2j>y9?3xSMC7jUM|0JpZGY9Cq;uNqcasw5=r~Sb)83=flov7D_1qX>I+yzJk zAq)<2KIsyK<|(9N_N^dz(X8HXE!2Ro5$cmRqn1W};F9pG0@Mo$fa;eYl*d;Qx$6c6 zc4g8apJ(Q)Q^*2mRvODc`G5bsz*jrp;!>HDm~oUD5Vqn+Sv^t!nO*mDs8*EWC~oEf zPD=~W^LJ6X_Yn6Vetg^hlJC{lg=z%)yC3c48}orSMN@y*6gL=RT})ms{dN7VvnRU*H8Tt%6L3VLVuKpLvGQ{23Wt)pgx{NdV*@uH75CCIM~z=C{Sr_c1<;czNDN zSs3I}Po93Zu#EUh7t9Tf%R@`tDfhnaKzLGTdQkZ(0={qhC>Eq;fIP?il8>}Jcvov5 z@2}tlG3tqv`~u2AJ#~P&Idy1${QJ01ME6Ip=Ts0aShu=; zuQ2%R%X+u0qYJVxCvICeY69_;a)E+AAIyf}UVN5g0aQhM-x3i2x(?p$>pQL6^z2~Q z$Y{?HxjwAi+u!+;(i4_G$OaDL20%lL7M}M?0H|>bf22R93B2wNSkz5FK92l2Gg~!( z1_N)hqI5%_1wq=~iT}bZBMq1gP#E&^~fdJJ&=yffCxpen^JI`X z*T{}Aj3G+QLs>IEa*(TZ>39>BDOmU{O(c+LgV7h($mBRZXu)CB+M_54Wu65gyL(t5 zQhrs5eozy72WwA{wTQxlFRKjaqHUm%@>5~w4k`$YEM;CoKL-{4O1C(V~_=qW%)Sjxp?l>yvx0|{1*e)~E= zp6@O>Hy=Kb$um7iwjc?$XN#G<_DBLDF@sf_ksJ_yA#y$AtpGyw#|N+FnSpB07k|ec ze87HutYFNkfAiGFa;7|2at)_2$jK}d^-P0i6TxSuG_&|hH#39Q-E zNsVrh;p*Xdvfm0;h?Jk0ezu446|Ka$9lx%37JMmc%Ij?O569EjakYhi3=?DR#Ivk(Lf+-X(Hnx zxiXk3j4U}mP51>LI;&O-Jq04GC#J~8J^qef3T)=2OE|9qvjQE z5Tr`beXWxY?h~G{R-O`p>+x8-Xx9Y5hqdm_**bcl@#3fEE5U(dBgBP)0}F_zR3LVG zzZpqB@#nTF>Va~S z8)cSOD^fJTR1xMJ4DPP%*)rZ5@aW?9hSZ0?V6^KRiJ+F~KgQQ@S>|dE+nYZ5`EPsO z*0@8|vi&nJQ1Og^$ayLPj2^+!63zT@zuiPMsZtO$z4$VtQ+dH`pkBR}P8MRBX=1ZX zB%m_xlf?;J~*;oZF_ge%JeAYHL8^5F-x&-nFBADlv9oSA{2Gm!}S_ zGA2}O`)$BbDD^?vNo(M6*kkf3IT+5CzGx}hX#}0@gxFf*OGu{em1$w^PspXV(ySI8 z9q^XPcWo3(0JHb^*~(m_fXYEZ$RIon0>8MgTZY8KU5Ru4H)(FZ$~1WO|iH$ z=ei4ohlv7RV-n%z-2ErDPL>cB^7v4SoI8|u-nYbZbq8J^-MNrwgs{iRK$1|<3$z~b zwmV*S0h`L<(Tua6|LgJgTjO!7{$Ck*7C_)~X@DPW_?m5&oLPa~?ZvHZ8DHS*vxpM# zaEFUuQkXo1@ktCTQ8NmxqaS&dVY1d$4qGBx^m zF(r}4VMn-v1QI*#xlLvc{n>BI&c9Opjt?T{_M`hrfBoI;M?7ukVEE{0uM3S~Ttcne zM1eUG8jM~t0)g(AUGK4H5RErvrn!Fd&|pGW`cfqVZsH%;owRa;XIjUr+CI1e;qEsK zxF#0xPNIAA7OosTCmY&NBxwfX^m~WK(cjl@9RrSw2VMf%~ycrxkCqU zJA*n-64d7!2~o%U!t2hrdM^h2fFotvDB^%OjJ=w;R&&Jy8Vsx$iyykcXMAh!ci;!p z+c^oy^ur;P@5K4W_;5IpeW2?F#|~IIwM_Kf$rJRQl9T;11L0jq?Gj0_0x&o>Wh{G1 z!vmbVOD`Q1KtQ-MFd={IcfTKvhkv)c_&no%lSy}oeIscYlwtu&<&O!@pLc;5sk&1g zFHAt0wyxOJUJEF(C$`6*R05~mB%*``1ma3Q_qTPiOK=d%yr zC$10yHOJzf^hN=&5@uCz^`(Y?*Uz7&eN~>+yaB5*=%ZdQF;oV?rCqPtk}0g9qV!hK z`|UR1N^P6^(A5wclWp&up0|KxHO3cey@s&1=WUD~?Thb@kKe}SSoRce<&NujeBER! z#u+Y`f~JM|koI7E;Bz=ZV4bQ8)xPESy*H&`wpk=CFu(vxWBGXlc>F;+hV39Fz6%6X zE6k-HvIDF9i}`aI%OeoIzIJy3$wWD+Jq*VX9K-MMf;V{F*vc&7aBs(#x&`!n zz2$n|BZrQ?FmR?nmxwqXD(VD3P2BeYw&XIBwnr|oJ6)k%+rb*W{vxjGa=taBJ=+&* zO=b)6?ly#@x9p%X@Z}knoz9RS4r^h@EaCi=SL9)N2MD^ex=dCX1b3ciHj16{feR)1 zBai8ApybWF2ir?+AoP@Yd(~?kaM2gL5?A&8?*im+`GTrfG~kI;sPzf=k4Q3)F&RUV z6|nW<`o!HchZ`AGB?4`BaPdxz>~YD#-{kvS{{LI#S0j@3j#u0Ougl4FjB7hIzOkMa zu5*HSNhivd&g()18EV>k?`It+YU&?Wfo-RT zhytoA9K600r zpli1+NU&;;zP<7N{`#Z17^n|Z94qqwjt8f=>o>aAe=Wb+uEx&!C=`sI-&8i=6AB0J z#}-!NIzqBPuO%04D3IpR9odI#4$-;AI)xe`FzlXJqsrk6=E^GC+X}_PE=FUo>P^J4`HKWjO^Re&4jePw006F@20l;2-A5Li_L54WKd(ASIG zUNn>kz%jCUpSB`TFxDt{h%*m^ywn}EXVL37-!SZ#p<~m9FP(tJC8Q0XmoFxstnz>t z=k`AvH(Ema5PbPk(~n4mHR<*KL%%-n+nCVU@w`G2g6YPZcb)zIx!$k5?>sGCQArW- z9S>jY|64J=JxC;bP7}^IpRRLWj)VHQZg7AWuRI*JaE@m7H-u}&5-i*@CUC&>;44?>fuH8Me2LA9*8A#jUl;#6Uf3f! z$IhYe|KSR5sXpay36J*~UejX=fS8Q^#}ZDse>dKen0f8VCQZI8*Q!70(w=7wwId(z z=eq;|>z)~%aJOHNBOkim&vaoM!K=XLi)s&p;1k(y^_LIapdvX>Xl^nbyA^xlyw9Sb5wKoG$KsGu|ph+sibK?DUwKm{A3 z2uKs9sE8<1Y;e~fbqpO2&pzjT=k9aA{r<_wo2(@3NoH0uGv||sO3`+bxa+2QV(^|R zdT#F2k8dIPT3LI2j!wBExa3g&^-f>>abwL;(eyf`3jLunE#gV|iI-oK|MH&H)a>(a zS(Ghomt``05$mq6{qd`Cfhzsz>kEJL_4yhe!Q!+`xR^x7-`u*d2XsRVdEb zrRR*eL}GpVx=lyY!$BHf{pE_`dLW0=k8>ml-tWrV=Lx>`2!%MCcdO*#qs;$gdx8j# zf7+3;R6`Nboi*#)h;#1b*oGRt6d6S2T@`EIB!s)QsT@zqIN|dk+g*A<0d}*J7kvo6 zj>OKkfQU=RzwN)T?f$EUM~1Q z>WDdDeBoz4%OCOJllC9bZH>dop~o4=^?tt3#;16)*3wiP6!L*HIqhOtyXyJMe5#d5 zc;OK(f0qjHGs7ZWgjoh{IMRv=nVEYLit!Qghz`St!AKJ}T|V zm*b2G8%ga|RZO6~*|y<01HlW`B5g-U@Sv)fwy32#v6RNqA4VSXGBcxej+x2c= zBrl$*-cYa>{oXqc;%FF@19HzTn{>bhQr)B&q2KBm;wh`W2rTYO^pQ8y0IO7vQDwa_ zBpYJg`YibH@^dYThw?lB;w1+U8DkqRG_7R1l)Z!xevCJpO+JVqas^ApjJhsDmFMjb zr@P|z-3#}(Ix>EflfQaD^cU!!P^GYi&{HjDs|!I$i+%m^e6l0v=;IkwQwV-8=k;-w zBE-79bk)LUYF|8?ve`*Xu?d}x3GE3P(b#cKda0dM5CRhV(pDY;d=nQEYN*&DKtXx1 z$;b^BzE9_`klgXOc;6oVLIVsWOYU&uIZLvww#u4`<^P6Lof+FxLHAX-m`>Tfev<;+ z2G|aMx*-R}K_{8(vm@W+K<`Re{Ia(rq_*hxv)iPBPj@&NFYysPscOMFmuETQY`3dr zzET}d_xBvXoJ;Vq$sbH?RS|%J1I0)TF>lkPd!F@NKSD}6Fc#UVz>9sC2UDw+2s`6^ zHl}JqFP!2s&a9-b3BJ-X^Bd1B;5yDvDYDHQe!K^_iUmcX{=kCK`qyi*j^bW>z;;iR z(2h*#7#m`=ijs##TpudCoW%VOYQfj1RNhS25~of3PrhzaL#xmv*W1&2xOsA0`sXEt zzLi2W{ZspA(oAbqb~241Y<&d>{oYyQX@0Krn3Nq#_Btrf7Fgf_N2PzcsMR-|{_5lD zKh?f_)s}OUU+l3zrTQ*$`?d4oSUa@X~ePzbtB({zXU!PPF4nc}n)f~=%i zoutfhAaF)w=`$@vksW%p=(+~4AWrBgZ^Okg@|=aV!|@n6%A*ZqZ0*Cuu0WT&xpC3y&Dbu<%5R_pM(`b3)~;b%0}ue zaesBM!9@kv<#PC9@t?NbIDKlN0%s2U_F_#AB2h4E}($;HtqVqVc2 zKVdtmg&-r>*7Lk_5GZ}aoZ4-GYWj_u+P6(H^3hYU{McG-n4}u?J3*}9OqR|xyjf4! z_gr;;hvl<__qu9kg6n=?QAZ*BsC4>~TzKcMcba zC$U~qpZY`F3)={tQGk&P6v^*?*&j^U6&sk-pGnrkbI~-?GEF&OJwCFQR+w-OI2kM+0=Nn^ZJ$G)*y2#1r(VR=p*b>h2^7new4 z_pmr`*@Or()Wk)ZPsh4HAyx`EFENFH3N$wABGttq!rUNgZml=5;R6nKvZ{fYSr{c!O5#oO7I`T^8 zCQaOa=osW`e2%2eYFj_e$A$u~&e5*hhRE0x&aqf#g;Qna8memYFp{d9IaVZzQSEaR z9VOCOTFn!DR9FzxXZq7PvKb(es&cUQg8;VmP|?iWiGZW*J+PR*#{j@S|voEdeY87c)vws;ND2+898 zJ?$O&ElglCeCJp0s(@PW(=|Opp4g~orupupAMB-Ot!7n%Ay;r~dRROHiPay~-w=Ag z6w5;+PK}{BQ1PJXNmLZ}H+*n@FzJP88)NN}O+IMlU~`W-?1DM#?bSNt4k)O&BGmfO z3heLJ-Z+tCi-U;3}x8Fben;*k&CZ|nO{jaQ0T zqm(Z5jWB<%a{DEPD=HQbaNJ$AgSm5_AzS1|=-9FpBnitF#;`(1&=y1*byhMn!B*qrwO5@{1 zF|k^BCHi^baBOKrsy_!v@t7?XTne7rApxsMesEk`qEKKe4`%BjK^IvgXle_ysjM)B z+mgj}&j1H-^p~1_8nXtEKwf(LK1*c3c6ADUZGw5}y>|=4N&oA(Uf9!;Q?gFzTRY@C z%=aAo@UwpJYyWci$2*dDH~jSdsaXEbYM@U9@&;#L<|Zn_^+ucu?{XEqtSZuurIf6-BS7L*u^f3Q&Lb*rETl6fA_+vZn2pN4=nw{}U=@ zxNQz;@+15OMAvQ9RCHFu8<}B+_58-D*(CZ_wQLP8m|s>+PjtjBo2+oFb-ajs>Fl;^ zt29Kl86H_ar9^0GN2uM%BB?G-uKdo46-b?DvoBMj22WR9Swt)?gs)3|(hF9GyPNyX z8Gk_1xJ+&0jR7(5EzTt6+ahXd#J-y|_NY6PC*4AvD+aGC8OmjILgs^{H!t5RW2LwP zd1je4de6@)HW!&=LhJUM#uh(NclSkaIphxVl74nQAs4LbDBeOI>4u8DF0R=c;`?s2 z$(%mX6047IXy{*MkJ}H>zA+i|D)mi^a~pgzaiEQ zm%1g-+5c?c{*RXPSIHM+EaNs56^8REX0qp}0`O^AX6_iwu-)w_-QvsjxShT37&pyY zP?0!=L=RgaAjN*%#!vy*Hb^o(=v0NCyuRO3q1Bivm{qy5(+JBeKk_^6=0om;j8lLi zGcr3SJE~6bLPOm4^w}G$K)dweoc%mK+#8l8B~r7XqGsn+MFRnR#e@7{Kptt81m+|I z3d#&sK+4$>7^p9f!?C&dDO%)kmtH~Cj&z2mB^L42R1ZJCa>Ecu2{ zvzS+BKCyo{psH$G@KPh}?bF-D)Hy)b74vv06*oRgg)nBbN}}@Yl$q03E^PX^!qaM< zGWHzOal9@fNSxnI=?DqPVIofC?s^Y#$gg(23 z)Va+@s9TfL_jtEHI&)XXI;!b_cdvX;Dsf-gaR!Ci6v75v>kn0^x%R#fGN zmz~Nzfsxa%EOlqP z*7V#E{T{;{4l0H?DYEC4(EBx@Qopm#gN8W2bPR?Pk$E465P5mAxbn%L=b8E)PXDg-JHf&IyBr+xnxs4V%;;x-h?9cT zGpRwYc*T4~omA|EcbnwYoktDfhoi0R%KBh)-cmhtRR+Im&($w=pUU(uk{*VyOrRn5 z=e9!1iCj)%$Vk@i-C@ju7|r|SZQPP*s-*Ov*d`1A@zCw>ITvIjW!E`yF-u4eC%=-ux&x#4ac7+LNGc!N7CXIC0v@xP#oN=_D zmEbvgWpU4+OF9^Ug!&M6Ax4p!=B87FE2F&PWxMJhv;j$hQoj$?Bl-xh{gTNN!6ZX&D5dZLUki{|CpVEKTZ~y53lwvh^ zQrM#m-=+chCEIVa(8NCaKN&Q5!0)p4IRa+>Xq7&b(!wsTj?fM$H! zg)LGMmCAH-$+NuZI1sD8z$*NwaaqoK=Y2J??+xwq=o6zhM(|y8hW0J`-^!Ub-*C0t z?5`{LU-AD(!_$l09Q`8ev0AIKe4fP<*`@nvNBkYJ_ly%$P5pO2_}O>M*q*(Rg@Uo4 z(KABtP99X_|A|-^RbF@HIPFT@d)vz`@xW*`s)QO79cDyelk-rgQC|s(g}mEk%0*$f z{BWi1VMUZIE-Snctbh=8(IM9WaWs>=_fiR{BVtm}wQZX`YJ?~{FLStHb@d6+%O-@M zQa;bj%4-6+@~}s((2H3A37Vvw*~b9!(}@dd)yqM7cA%%~-b#oL-tBbmqC;x|-{hQA9qG21be1{! z=(sViHXF{2Cg@>rQa}M|Obdav^hFiY#652mymTDzR1qs_#Q4QU9ZPjLnjV!_!HbS4 z+4ufZ(9nM@J52ESi%XXnjJ6VX%c6aadKcAT8Tg)id8!y%_uXHT>!pZ05`7ZKxxJ9F zkW@LoJrKrym4StQB`Z8`4nI3u?q$eT2S@5zRdzzuX0yd#cM-U#(%Fo+eY+^a`V}91mLmM z`#<)E673!LWeiRZLRGCNr{)n)z+p3Hc0y7(YB zK(M*kA{<&r=@UPYtpmqwzE_+fAS?3VaWS(rMnXu*`}6p~3x0<(f|vVo@SSuw8Ub*m zf65+BW8%t^`7=xSb_NZp0K#42R%P#NkeK{Jim9wIO+ut zE*;-VaW|RRXL}*?lD`6CO$Ro+)qeM*2<_o*T9l;$*X**U`{mMzP4~+ZT_K25#i|F% zvm|h@UdYJ*k{sAJ<_cVQP{z4?t;dwJq|x^hWg>9xvWZPHWM zxhw=tzhtsgMSKWl6kO;U=7GmXE8gal+!(ekt$FvE1FYq?Vab}JzVvanp zPWr-yTfuuy)z~B@zNrn6n9^*X9vJ9ri%Sj+a(aT(>xJtZQD@kG3ah7=wFjpjC&g{z{~ogXvAe;+D1361P4s#s zQhWB$y|7si1Ims>eu5Y2dHo$zb}~!Iz4k~yL#zwvFPtjq(04&r=E|t8m7Yj{T0Xt- z+!H6w_8$L(+7%*dbuRgrop6nv|Db4v4JOj!2gL%d2)n198x*+oVRz52uhU`;MC`F* z;+YA=+at&O#O+WkoRHk(Y5i?HMxLh6eYeaS#rr8}pGMeW($|GPDbN~%C;L3z*Bikr z%4N}&(gpqdOxjuIy-+vxg!;)=!XLapB(I9lw;t)Ry}H2a3B}PP3|5a_u;mE(<;fOT zJkeJ2_a*Gm?N3QE*|odkaqB{o9=QWT-p459!5g$|^+EaeivM40C zHicv!Rm2qKapt2<1aF4q8roWopY2$_%I__WDazJkYw&W(2$l2$Z_phmQP|VU3xfw) zefjPb7_ax(f@&Pz4>d?A9(85x9sRZwzh;*Lgjn4jEn z+G5TGXV{kA+8S?*eRNvcT?b4EJK%Kr_8cX!DL=Y5Zp{9zUC*-`4Lek|zs1{v9VClm z_%R^((D`-`p-+s+33?^527|RqOw3Q@u=1_nY+An>cE>7feTh}a0tM+q=W$&adMzJK z-@=cY+i`qdee_VKa#&HTBLwGL*Xffhxe=E8Uf}8jL7Xvo$M#HWHQ2Y0Uw_Z4hW1q? z$PD|SHrGsc;dTHvOAY5{4|^l6z526koF{^xZ{l|&{7Z;=dh6*lZS0$#pK@|nh5u(M znTSjS=!$>&z|U?7|6?vg0{#Zb7~bh-xLXgF+kDanh;h0|yZc6Gkuuad&gW(mIs9gc*>gUWw%g{jvlIR* zbtN5%(b}+H)LK{S&Vo2U&wQN;dgvIbrP?l+#*19$Yp0lHuwJ8RNv@(b4sPFXH=5vu z)7uV8ly$g+=EUPo)p9Cmwr8+Ew8jm`M3Z2Hp;zkKDtBitejL`rW3kV9zS+8@I|!h$9 zg*DcG_)yk#&)(P*<_Sj_Ek^Con!nI%ZtZ}>w(nKeb*(|erDfexi`tOkI=A06TL-oG z8miv=s$&^5#TZMVE-J}}?~rrpAcE!E!uEC@n0lu>9?npMsIN{#Wg-nGSoWmqpI|`Q zCSNv%jGy_t2B;+GH;Yq1{-x>_X@^!4?;EnnvIRBpG}M<3&WJ)c;;HxinI4k*{=ipa zoQ9y0I(nDZ{r$J^|C>QpjuU#cST$ExB%YXZl z75vf|I8hz3wz0wlt53vag*fWKu~Kj%v%nC!58B%jUJ!ift1nuykISR8kVe``MILd+ z77V2ul`$t?zHfK6AtdAIXK01=(Y&he&HNWLIPPj>+c?b)+07;VZc9Z_zBht*F^LhH z4KC|n*NVZ_x@MBRni2lTC+Zia z_HF+B5y#)fFKo=(mYMn8zVmnae@y2YzR(jzeU)UQ`dRb5R?|0m`7z%=8{a=!C7mK; zf-S0vrG+csk?P{#mU3(mfF7Bh4xgwkGIgd(-tJXJl3qv$1)~zOn=;0oJmiq)=x?kp z{at@;A64HUNt|yF=UC(&A^d}OKKJZaQ~a(!Y0P5ai{GG*ppAhM#;ez$d2&q^MTR!= z26ULU6s!>18mJ}}?*yhIsV~%m;doeE%vE6UeSY}e`yI$-Y?Qpgfzd~a`BrnDF!-1! zGT|qPu6t`=%O0UaW7GM!Z8b~4#9G|gEAWBzQsC5Tb8-s!Ns65_+DHYxbm^>#C7Kv8 z+Fig}9*r+yO>N_CDk+bh{UBB?H6u)?VB7F*qo`nnN>#m zV@>3B&os|(|Np8#s-7s*el(*5L)_@IYQBrV&Ht-%{;&Ms*DScT!CD5%Q(}FXGr-L> zd#V-n1_)c;9y@f29>06fd8S0^?Lx~yALGg3)Cy`>G8CuqEYNV1W{^|Xk62FD7l**@-2^E z$4e+7RG0TNhoKi5T;9D=Ke7>Xqg9nDrjh@2`~K`}jWiY38+m^B)A_T{|GV zbbTa|J-9$8((tpM`QKGf{^KQF=jkOOw);?*&23}!`lg;-b4nAgC#g=0&`W?}bCgAI z$<_bw<@Dc$_pizczm>>VzIPdp_V23Qs}1~^k1uAg@ioa0KlKAczJN4;^)$)gaB*SE zAA}y+PPL(L&J3GveJ&)Z5qwqNYwp*iS)wx7a^rr2*Tv$vxo)Ji52%;YG}ROP=`n^; z5mnQl&sCQRCUoCUm&Jvf+)N^b-B=5c*oD4edtB3F7HYD!z!p!3%&f@OIFdFS%wQ)3 z125eJI#l#9;=j4ZvWXmf2QmT%*%wJ9Q;JK>bxW|SZCBHsqG?jym1yctu4z(NZAD9V zFE!Gf6*Xjzup!#?ty0)YKC~4)=6bh80v!z*UDDkicr)9ScyzZD&iLLF6nm)q6D~jE z^k?JtC+qsuygYD5dk5Di3oj_8xSrw9@k4~jKKWPMR1u_LJlYrO3GJST$=UO!FwuK1 zQ|kC@z3HE=|38&qQn>fwR?#n{-Y&|GOdDL#ZAxu>zv92&p5{M=-(S6+omZ!MuZ(^m zy{GhNmN}^nD$8{_&a8o$3Z1YuqY6Op^1h22Z+x&fw(6Wwx--G|Oh;z&EC`YHnXed& z!cnD|b0YddFq|4rQOPC!e9m`kr!JM^DsxaWz3d%%OUV6dEr+IXDa6!X^L{+bg5}hw zR93sM1g%g=Jq;ftr086?g|A^iVoCM(g{0cu^9Vd3M1TM&Ymc+F}FIcs09&)nNmsY#zG>ZwDIaEmVggEjTMUfg%{gtji2|l6>$KFnDUW z4z=-XA{EvKV1Jt3!_3C@SngR6xai>k4GsFv src/shared/topography/topography.F90) +expects *raw high-resolution elevation* and computes subgrid standard deviation +itself by binning many fine-grid points into each destination grid cell. ERA5's +sdor is already a precomputed standard-deviation field (ECMWF derives it the same +way, from a much finer DEM), so feeding it through 'computed' would compute the +"stdev of a stdev field" -- meaningless. Instead this script regrids sdor +directly onto the model grid and mg_drag reads it verbatim via +source_of_sgsmtn='input' (mg_drag_init's INPUT/mg_drag.res.nc restart-style +read, variable name 'ghprime'). + +Regridding: bilinear interpolation onto Isca's Gaussian T42 latitudes (via +numpy.polynomial.legendre.leggauss). Interpolation is done by hand with two +separable 1D numpy.interp passes (longitude then latitude) rather than +xarray's/scipy's interp -- this works cleanly here because both the ERA5 +source grid and the T42 target grid are regular in longitude and monotonic in +latitude, so a separable linear interpolation is exact for bilinear +interpolation on a rectilinear grid. Spherical-harmonic filtering is +deliberately not used: sdor is a standard-deviation field and must stay >= 0 +everywhere, and spectral truncation can ring negative near sharp gradients +(e.g. the Himalayas); linear interpolation of a non-negative field stays +non-negative. + +Latitude convention: mg_drag_init reads this file via FMS's read_data(), which +loads the raw array by index position into the model's internal field with no +coordinate-based matching -- it just trusts the array is already ordered to +match the model's own domain-decomposed grid, exactly like every other Isca +input file (era-spectral7_T42_64x128.out.nc is ascending, -90 -> 90). Writing +descending here would silently mirror the whole field north-south when +mg_drag reads it at runtime, so gaussian_grid() below returns ascending +latitudes to match. +""" +import os + +import numpy as np +import xarray as xr +from numpy.polynomial.legendre import leggauss + +from isca import GFDL_BASE + +ERA5_SDOR_PATH = '/home/links/sit204/gwd_tests/69eb30085e940a6d80b1ce305b746b0b.nc' +OUT_PATH = os.path.join(GFDL_BASE, 'exp/test_cases/mg_drag/input/mg_drag.res.nc') + +NLAT_T42 = 64 +NLON_T42 = 128 + + +def gaussian_grid(n_lat: int, n_lon: int) -> tuple[np.ndarray, np.ndarray]: + """Ascending (-90 -> 90) Gaussian latitudes, matching Isca's own grid + convention -- see the 'Latitude convention' note above for why this matters.""" + gauss_points, _ = leggauss(n_lat) + lats = np.arcsin(gauss_points) * (180.0 / np.pi) + lons = np.linspace(0, 360, n_lon, endpoint=False) + return lats, lons + + +def bilinear_regrid(data, lat_src, lon_src, lat_tgt, lon_tgt): + """Separable linear regrid (no scipy). lat_src must be ascending; + lon_src must be ascending and span >= [lon_tgt.min(), lon_tgt.max()].""" + # interpolate along longitude for every source latitude row + by_lon = np.empty((data.shape[0], len(lon_tgt))) + for i in range(data.shape[0]): + by_lon[i, :] = np.interp(lon_tgt, lon_src, data[i, :]) + # then interpolate along latitude for every target longitude column + out = np.empty((len(lat_tgt), len(lon_tgt))) + for j in range(len(lon_tgt)): + out[:, j] = np.interp(lat_tgt, lat_src, by_lon[:, j]) + return out + + +def main(): + era5 = xr.open_dataset(ERA5_SDOR_PATH).squeeze() + sdor = era5['sdor'].sortby('latitude') # ERA5 ships lat descending; need ascending for np.interp + + lat_src = sdor['latitude'].values + lon_src = sdor['longitude'].values + data_src = sdor.values + + target_lat, target_lon = gaussian_grid(NLAT_T42, NLON_T42) + regridded = bilinear_regrid(data_src, lat_src, lon_src, target_lat, target_lon) + + assert regridded.shape == (NLAT_T42, NLON_T42) + assert regridded.min() >= 0.0, 'regridded sdor went negative -- unexpected for linear interp of a >=0 field' + + ghprime = xr.DataArray( + regridded.astype('float64'), + dims=('lat', 'lon'), + coords={'lat': target_lat, 'lon': target_lon}, + name='ghprime', + attrs={'long_name': 'subgrid-scale orography standard deviation (from ERA5 sdor)', 'units': 'm'}, + ) + ghprime.to_netcdf(OUT_PATH) + print(f'Wrote {OUT_PATH}: ghprime min={float(ghprime.min()):.1f} max={float(ghprime.max()):.1f} ' + f'mean={float(ghprime.mean()):.2f} m, shape={ghprime.shape}') + + +if __name__ == '__main__': + main() From cfb5983c6e136b872cd61f05764271144aae1610 Mon Sep 17 00:00:00 2001 From: sit23 Date: Thu, 13 Aug 2026 16:09:23 +0100 Subject: [PATCH 4/5] Fix column model build: mg_drag.f90 -> mg_drag.F90, guard grid_domain import CI failed with "Cannot open module file 'transforms_mod.mod'" when compiling mg_drag for the column model. Root cause: mkmf's generated Makefile only passes $(CPPDEFS) (the -D macros, including -DCOLUMN_MODEL) to .F90 (uppercase) source files -- lowercase .f90 files are compiled without any -D flags at all, regardless of which codebase is being built. The previous commit's #ifdef COLUMN_MODEL guard around mg_drag.f90's grid_domain import could therefore never see COLUMN_MODEL defined, always fell through to the transforms_mod branch, and transforms.F90 isn't part of the column model's path_names (it uses spec_mpp_mod instead, which provides the same grid_domain type) -- hence the missing-module error. Renamed mg_drag.f90 -> mg_drag.F90 (matching hs_forcing.F90's existing COLUMN_MODEL/spec_mpp_mod-vs-transforms_mod pattern, which works because that file is already uppercase) and updated the six path_names files that reference it (isca, socrates, socrates_column, column, dry, grey). Verified: IscaCodeBase (RRTM), SocratesCodeBase, and ColumnCodeBase all compile cleanly. Re-ran the mg_drag restart-domain-fix smoke test after the rename to confirm behaviour is unchanged. --- src/atmos_param/mg_drag/{mg_drag.f90 => mg_drag.F90} | 6 +++++- src/extra/model/column/path_names | 2 +- src/extra/model/dry/path_names | 2 +- src/extra/model/grey/path_names | 2 +- src/extra/model/isca/path_names | 2 +- src/extra/model/socrates/path_names | 2 +- src/extra/model/socrates_column/path_names | 2 +- src/extra/python/scripts/regrid_era5_sdor_to_t42.py | 2 +- 8 files changed, 12 insertions(+), 8 deletions(-) rename src/atmos_param/mg_drag/{mg_drag.f90 => mg_drag.F90} (99%) diff --git a/src/atmos_param/mg_drag/mg_drag.f90 b/src/atmos_param/mg_drag/mg_drag.F90 similarity index 99% rename from src/atmos_param/mg_drag/mg_drag.f90 rename to src/atmos_param/mg_drag/mg_drag.F90 index 99f3fbef1..a6bd852a4 100644 --- a/src/atmos_param/mg_drag/mg_drag.f90 +++ b/src/atmos_param/mg_drag/mg_drag.F90 @@ -17,7 +17,11 @@ module mg_drag_mod set_domain use fms_io_mod, only: get_restart_io_mode use constants_mod, only: Grav, Kappa, RDgas, cp_air - use transforms_mod, only: grid_domain +#ifdef COLUMN_MODEL + use spec_mpp_mod, only: grid_domain +#else + use transforms_mod, only: grid_domain +#endif !----------------------------------------------------------------------- implicit none diff --git a/src/extra/model/column/path_names b/src/extra/model/column/path_names index 6912d90f4..62758fe8b 100644 --- a/src/extra/model/column/path_names +++ b/src/extra/model/column/path_names @@ -105,7 +105,7 @@ atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_taumol.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_vrtqdr.f90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/dry/path_names b/src/extra/model/dry/path_names index f7caddc6c..890ac6372 100644 --- a/src/extra/model/dry/path_names +++ b/src/extra/model/dry/path_names @@ -25,7 +25,7 @@ atmos_param/frierson_monin_obukhov/frierson_monin_obukhov.F90 atmos_param/dry_convection/dry_convection.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/grey/path_names b/src/extra/model/grey/path_names index 8edec5d9f..d5a045000 100644 --- a/src/extra/model/grey/path_names +++ b/src/extra/model/grey/path_names @@ -26,7 +26,7 @@ atmos_param/dry_convection/dry_convection.f90 atmos_param/ras/ras.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/isca/path_names b/src/extra/model/isca/path_names index abcba6b80..74046e3f8 100644 --- a/src/extra/model/isca/path_names +++ b/src/extra/model/isca/path_names @@ -101,7 +101,7 @@ atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_taumol.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_vrtqdr.f90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/socrates/path_names b/src/extra/model/socrates/path_names index cab81e74d..f7762fb85 100644 --- a/src/extra/model/socrates/path_names +++ b/src/extra/model/socrates/path_names @@ -26,7 +26,7 @@ atmos_param/frierson_monin_obukhov/frierson_monin_obukhov.F90 atmos_param/dry_convection/dry_convection.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/socrates_column/path_names b/src/extra/model/socrates_column/path_names index 20fae2c23..f3d3f4ef7 100644 --- a/src/extra/model/socrates_column/path_names +++ b/src/extra/model/socrates_column/path_names @@ -30,7 +30,7 @@ atmos_param/frierson_monin_obukhov/frierson_monin_obukhov.F90 atmos_param/dry_convection/dry_convection.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/python/scripts/regrid_era5_sdor_to_t42.py b/src/extra/python/scripts/regrid_era5_sdor_to_t42.py index 3219bbdf5..a4e0890f9 100644 --- a/src/extra/python/scripts/regrid_era5_sdor_to_t42.py +++ b/src/extra/python/scripts/regrid_era5_sdor_to_t42.py @@ -12,7 +12,7 @@ regeneratable from CDS) -- edit ERA5_SDOR_PATH below to point at your own copy. Why 'input' mode instead of 'computed': mg_drag_nml's source_of_sgsmtn='computed' -path (src/atmos_param/mg_drag/mg_drag.f90 -> src/shared/topography/topography.F90) +path (src/atmos_param/mg_drag/mg_drag.F90 -> src/shared/topography/topography.F90) expects *raw high-resolution elevation* and computes subgrid standard deviation itself by binning many fine-grid points into each destination grid cell. ERA5's sdor is already a precomputed standard-deviation field (ECMWF derives it the same From 773a532e8e703913d1d18baaca05fd7255c2b990 Mon Sep 17 00:00:00 2001 From: sit23 Date: Thu, 13 Aug 2026 17:27:11 +0100 Subject: [PATCH 5/5] Fix trip test comparison crashing on column model's degenerate lonb coordinate conduct_comparison_on_test_case compared two commits' output by subtracting one xarray Dataset from the other, which aligns both datasets by coordinate value before diffing. The column model's 'lonb' coordinate is two identical near-zero placeholder values (a 1-column grid has no real longitude extent), and pandas' Index requires unique values to align/reindex on -- so the subtraction raised ValueError: cannot reindex or align along dimension 'lonb' because the (pandas) index has duplicate values, for column_test specifically (every other registered test case has well-behaved coordinates and never hit this). Both datasets being compared are the same test case run at the same resolution by construction (only the commit differs), so they already share the same grid -- a positional value comparison is exactly what's wanted, and doesn't need coordinate alignment at all. Switched to comparing each variable's raw .values arrays directly, which sidesteps xarray's alignment machinery entirely. Verified against already-completed test data: column_test now compares cleanly (no crash, correctly reports pass -- f48bfa1b and cfb5983c give identical output), and held_suarez's result is unchanged (still passes), confirming the new comparison is equivalent for well-behaved coordinates. --- .../trip_test/trip_test_functions.py | 20 +++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/exp/test_cases/trip_test/trip_test_functions.py b/exp/test_cases/trip_test/trip_test_functions.py index 53223fdff..f7ecae1c0 100644 --- a/exp/test_cases/trip_test/trip_test_functions.py +++ b/exp/test_cases/trip_test/trip_test_functions.py @@ -372,13 +372,21 @@ def conduct_comparison_on_test_case(base_commit, later_commit, test_case_name, r base_commit_dataset = xar.open_dataset(data_dir_dict[base_commit] +'/run0001/'+diag_file_entry+'.nc', decode_times=False) later_commit_dataset = xar.open_dataset(data_dir_dict[later_commit]+'/run0001/'+diag_file_entry+'.nc', decode_times=False) - diff = later_commit_dataset - base_commit_dataset - - #Check each of the output variables for differences - for var in diff.data_vars.keys(): - maxval = np.abs(diff[var]).max() + #Compare raw array values position-by-position rather than taking an + #xarray Dataset difference. Dataset subtraction aligns the two datasets + #by coordinate value first, which raises if a coordinate has duplicate + #values -- e.g. the single-column model's 'lonb' is two identical + #(near-zero placeholder) values, since a 1-column grid has no real + #longitude extent. Both datasets are the same test case run at the same + #resolution by construction (only the commit differs), so they already + #share the same grid and a positional comparison is exactly what's + #wanted here. + for var in base_commit_dataset.data_vars.keys(): + base_vals = base_commit_dataset[var].values + later_vals = later_commit_dataset[var].values + maxval = np.abs(later_vals - base_vals).max() if maxval !=0.: - print('Test failed for '+var+' max diff value = '+str(maxval.values)) + print('Test failed for '+var+' max diff value = '+str(maxval)) test_pass = False base_experiment_input_nml = f90nml.read(data_dir_dict[base_commit] +'/run0001/input.nml')