module cg_drag_mod use mpp_mod, only: input_nml_file use fms_mod, only: fms_init, mpp_pe, mpp_root_pe, & file_exist, check_nml_error, & error_mesg, FATAL, WARNING, NOTE, & close_file, open_namelist_file, & stdlog, write_version_number, & read_data, write_data use fms_io_mod, only: register_restart_field, restart_file_type use fms_io_mod, only: save_restart, restore_state, get_mosaic_tile_file use time_manager_mod, only: time_manager_init, time_type use diag_manager_mod, only: diag_manager_init, & register_diag_field, send_data use constants_mod, only: constants_init, PI, RDGAS, GRAV, CP_AIR, & SECONDS_PER_DAY #ifdef COL_DIAG use column_diagnostics_mod, only: column_diagnostics_init, & initialize_diagnostic_columns, & column_diagnostics_header, & close_column_diagnostics_units #endif !------------------------------------------------------------------- implicit none 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]. !-------------------------------------------------------------------- !--------------------------------------------------------------------- !----------- ****** VERSION NUMBER ******* --------------------------- character(len=128) :: version = '$Id: cg_drag.F90,v 20.0 2013/12/13 23:09:00 fms Exp $' character(len=128) :: tagname = '$Name: tikal $' !--------------------------------------------------------------------- !------- interfaces -------- public cg_drag_init, cg_drag_calc, cg_drag_end, cg_drag_restart, & cg_drag_time_vary, cg_drag_endts private read_nc_restart_file, gwfc !--- for netcdf restart type(restart_file_type), pointer, save :: Cg_restart => NULL() type(restart_file_type), pointer, save :: Til_restart => NULL() logical :: in_different_file = .false. integer :: vers, old_time_step !wfc++ Addition for regular use integer, allocatable, dimension(:,:) :: source_level real, allocatable, dimension(:,:) :: source_amp real, allocatable, dimension(:,:,:) :: gwd_u, gwd_v !wfc-- !-------------------------------------------------------------------- !---- namelist ----- integer :: cg_drag_freq=0 ! calculation frequency [ s ] integer :: cg_drag_offset=0 ! offset of calculation from 00Z [ s ] ! only has use if restarts are written ! at 00Z and calculations are not done ! every time step 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 ! [ Pa ] 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 ! [ m/s ] ! previous values: 0.6 real :: Bt_0=.003 ! sum across the wave spectrum of ! the magnitude of momentum flux, ! divided by density [ m^2/s^2 ] real :: Bt_aug=.000 ! magnitude of momentum flux divided by density real :: Bt_nh=.003 ! magnitude of momentum flux divided by density (SH limit ) real :: Bt_sh=.003 ! magnitude of momentum flux divided by density (SH limit ) real :: Bt_eq=.000 ! magnitude of momentum flux divided by density (equator) real :: Bt_eq_width=4.0 ! scaling for width of equtorial momentum flux (equator) real :: phi0n = 30., phi0s = -30., dphin = 5., dphis = -5. 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 ! 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 ! columns [degrees, -90. -> 90. ] 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, & 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 !-------------------------------------------------------------------- !-------- public data ----- !-------------------------------------------------------------------- !------ private data ------ !-------------------------------------------------------------------- ! list of restart versions readable by this module. !-------------------------------------------------------------------- integer, dimension(3) :: restart_versions = (/ 1, 2, 3 /) ! v1 : ! 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 ! [ m/s^2 ] ! ked effective eddy diffusion coefficient resulting from ! gravity wave forcing [ m^2/s ] ! !-------------------------------------------------------------------- !wfc++ not needed if calcucate_ked is removed. !!!!rjw real, dimension(:,:,:), allocatable :: gwd, ked !wfc-- !-------------------------------------------------------------------- ! these are the arrays which define the gravity wave source spectrum: ! ! c0 gravity wave phase speeds [ m/s ] ! kwv horizontal wavenumbers of gravity waves [ /m ] ! k2 squares of wavenumbers [ /(m^2) ] ! !------------------------------------------------------------------- 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 ! k index of the gravity wave source level at ! the equator in a standard atmosphere !--------------------------------------------------------------------- ! variables which control module calculations: ! ! cgdrag_alarm time remaining until next cg_drag calculation [ s ] ! !--------------------------------------------------------------------- integer :: cgdrag_alarm !--------------------------------------------------------------------- ! variables used with column diagnostics: ! ! diag_units output unit numbers ! 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 ? ! 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 logical :: column_diagnostics_desired=.false. 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 !--------------------------------------------------------------------- ! variables for netcdf diagnostic fields. !--------------------------------------------------------------------- integer :: id_kedx_cgwd, id_kedy_cgwd, id_bf_cgwd, & id_gwfx_cgwd, id_gwfy_cgwd real :: missing_value = -999. character(len=7) :: mod_name = 'cg_drag' logical :: module_is_initialized=.false. !------------------------------------------------------------------- !------------------------------------------------------------------- contains !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! PUBLIC SUBROUTINES ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !#################################################################### subroutine cg_drag_init (lonb, latb, pref, Time, axes) !------------------------------------------------------------------- ! cg_drag_init is the constructor for cg_drag_mod. !------------------------------------------------------------------- !------------------------------------------------------------------- real, dimension(:,:), intent(in) :: lonb, latb real, dimension(:), intent(in) :: pref integer, dimension(4), intent(in) :: axes type(time_type), intent(in) :: Time !------------------------------------------------------------------- !------------------------------------------------------------------- ! intent(in) variables: ! ! lonb 2d array of model longitudes on cell corners [radians] ! latb 2d array of model latitudes at cell corners [radians] ! pref array of reference pressures at full levels (plus ! surface value at nlev+1), based on 1013.25hPa pstar ! [ Pa ] ! Time current time (time_type) ! axes data axes for diagnostics ! !------------------------------------------------------------------ !------------------------------------------------------------------- ! local variables: integer :: unit, ierr, io, logunit integer :: n, i, j, k integer :: idf, jdf, kmax real :: pif = 3.14159265358979/180. ! real :: pif = PI/180. ! real, allocatable :: lat(:,:) real :: lat(size(lonb,1) - 1, size(latb,2) - 1) !------------------------------------------------------------------- ! 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 ! jdf number of j points on this processor ! kmax number of k points on this processor ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! if routine has already been executed, return. !--------------------------------------------------------------------- if (module_is_initialized) return !--------------------------------------------------------------------- ! verify that all modules used by this module have been initialized. !--------------------------------------------------------------------- call fms_init call time_manager_init call diag_manager_init call constants_init #ifdef COL_DIAG call column_diagnostics_init #endif SKIP !--------------------------------------------------------------------- ! read namelist. !--------------------------------------------------------------------- #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=cg_drag_nml, iostat=io) ierr = check_nml_error(io,"cg_drag_nml") #else if (file_exist('input.nml')) then unit = open_namelist_file ( ) ierr=1; do while (ierr /= 0) read (unit, nml=cg_drag_nml, iostat=io, end=10) ierr = check_nml_error (io, 'cg_drag_nml') enddo 10 call close_file (unit) endif #endif !--------------------------------------------------------------------- ! write version number and namelist to logfile. !--------------------------------------------------------------------- call write_version_number(version, tagname) logunit = stdlog() 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 ! domain on this processor, kmax is the number of model layers. !------------------------------------------------------------------- kmax = size(pref(:)) - 1 jdf = size(latb,2) - 1 idf = size(lonb,1) - 1 allocate( source_level(idf,jdf) ) allocate( source_amp(idf,jdf) ) ! allocate( lat(idf,jdf) ) !-------------------------------------------------------------------- ! define the k level which will serve as source level for the grav- ! ity waves. it is that model level just below the pressure specif- ! ied as the source location via namelist input. !-------------------------------------------------------------------- do k=1,kmax if (pref(k) > source_level_pressure) then klevel_of_source = k exit endif end do do j=1,jdf do i=1,idf lat(i,j)= 0.5*( latb(i,j+1)+latb(i,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)) end do end do source_level = MIN (source_level, kmax-1) ! deallocate( lat ) !--------------------------------------------------------------------- ! determine if column diagnostics are desired from this module. if ! so, set a flag to so indicate. !--------------------------------------------------------------------- num_diag_pts = num_diag_pts_ij + num_diag_pts_latlon if (num_diag_pts > 0) then column_diagnostics_desired = .true. endif !--------------------------------------------------------------------- ! if column diagnostics are desired, check that array dimensions are ! sufficiently large for the number of requests. !--------------------------------------------------------------------- #ifdef COL_DIAG if (column_diagnostics_desired) then if (num_diag_pts > MAX_PTS) then call error_mesg ( 'cg_drag_mod', & ' must reset MAX_PTS or reduce number of diagnostic points', & FATAL) endif !--------------------------------------------------------------------- ! allocate arrays needed for column diagnostics. !--------------------------------------------------------------------- allocate (do_column_diagnostics (jdf) ) allocate (diag_units (num_diag_pts) ) allocate (diag_lon (num_diag_pts) ) allocate (diag_lat (num_diag_pts) ) allocate (diag_i (num_diag_pts) ) 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 ! 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, & lon_coords_gl, lonb(:,1), latb(1,:), do_column_diagnostics, & diag_lon, diag_lat, diag_i, diag_j, diag_units) endif #endif !--------------------------------------------------------------------- ! define the number of waves in the gravity wave spectrum, and define ! an array of their speeds. They are defined symmetrically around ! c = 0.0 m/s. !--------------------------------------------------------------------- nc = 2.0*cmax/dc + 1 allocate ( c0(nc) ) 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. !-------------------------------------------------------------------- allocate ( kwv(nk) ) allocate ( k2 (nk) ) do n=1,nk kwv(n) = 2.*PI/((30.*(10.**n))*1.e3) k2(n) = kwv(n)*kwv(n) end do !-------------------------------------------------------------------- ! initialize netcdf diagnostic fields. !------------------------------------------------------------------- id_bf_cgwd = & register_diag_field (mod_name, 'bf_cgwd', axes(1:3), Time, & 'buoyancy frequency from cg_drag', ' /s', & missing_value=missing_value) id_gwfx_cgwd = & register_diag_field (mod_name, 'gwfx_cgwd', axes(1:3), Time, & 'gravity wave forcing on mean flow', & 'm/s^2', missing_value=missing_value) id_gwfy_cgwd = & register_diag_field (mod_name, 'gwfy_cgwd', axes(1:3), Time, & 'gravity wave forcing on mean flow', & 'm/s^2', missing_value=missing_value) id_kedx_cgwd = & register_diag_field (mod_name, 'kedx_cgwd', axes(1:3), Time, & 'effective eddy viscosity from cg_drag', 'm^2/s', & missing_value=missing_value) id_kedy_cgwd = & register_diag_field (mod_name, 'kedy_cgwd', axes(1:3), Time, & 'effective eddy viscosity from cg_drag', 'm^2/s', & missing_value=missing_value) !-------------------------------------------------------------------- ! 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) ) !-------------------------------------------------------------------- ! if present, read the restart data file. !--------------------------------------------------------------------- if (size(restart_versions(:)) .gt. 2 ) then call cg_drag_register_restart endif if (file_exist('INPUT/cg_drag.res.nc')) then call read_nc_restart_file elseif (file_exist('INPUT/cg_drag.res')) then call error_mesg ( 'cg_drag_mod', 'Native restart capability has been removed.', & FATAL) !------------------------------------------------------------------- ! if no restart file is present, initialize the gwd field to zero. ! define the time remaining until the next cg_drag calculation from ! the namelist inputs. !------------------------------------------------------------------- else gwd_u(:,:,:) = 0.0 gwd_v(:,:,:) = 0.0 if (cg_drag_offset > 0) then cgdrag_alarm = cg_drag_offset else cgdrag_alarm = cg_drag_freq endif endif vers = restart_versions(size(restart_versions(:))) old_time_step = cgdrag_alarm !--------------------------------------------------------------------- ! mark the module as initialized. !--------------------------------------------------------------------- module_is_initialized = .true. !--------------------------------------------------------------------- end subroutine cg_drag_init !#################################################################### subroutine cg_drag_time_vary (delt) real , intent(in) :: delt !--------------------------------------------------------------------- ! decrement the time remaining until the next cg_drag calculation. !--------------------------------------------------------------------- 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 ! forcing. !-------------------------------------------------------------------- if (cgdrag_alarm <= 0 ) then cgdrag_alarm = cgdrag_alarm + cg_drag_freq endif 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 ! the desired output fields, and saves the values for later retrieval ! if they are not calculated on every timestep. ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- integer, intent(in) :: is, js real, dimension(:,:), intent(in) :: lat real, dimension(:,:,:), intent(in) :: pfull, zfull, temp, uuu, vvv type(time_type), intent(in) :: Time real , intent(in) :: delt real, dimension(:,:,:), intent(out) :: gwfcng_x, gwfcng_y !------------------------------------------------------------------- ! intent(in) variables: ! ! 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 ] ! zfull height at model full levels [ m ] ! temp temperature at model levels [ deg K ] ! uuu zonal wind [ m/s ] ! vvv meridional wind [ m/s ] ! Time current time, needed for diagnostics [ time_type ] ! delt physics time step [ s ] ! ! intent(out) variables: ! ! gwfcng_x time tendency for u eqn due to gravity-wave forcing ! [ m/s^2 ] ! gwfcng_y time tendency for v eqn due to gravity-wave forcing ! [ m/s^2 ] ! !------------------------------------------------------------------- !------------------------------------------------------------------- ! local variables: real, dimension (size(uuu,1), size(uuu,2), size(uuu,3)) :: & dtdz, ked_gwfc_x, ked_gwfc_y real, dimension (size(uuu,1),size(uuu,2), 0:size(uuu,3)) :: & zzchm, zu, zv, zden, zbf, & gwd_xtnd, ked_xtnd, & gwd_ytnd, ked_ytnd integer :: iz0 logical :: used real :: bflim = 2.5E-5 integer :: ie, je integer :: imax, jmax, kmax integer :: i, j, k !------------------------------------------------------------------- ! local variables: ! ! dtdz temperature lapse rate [ deg K/m ] ! 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 ! [ m/s^2 ] ! ked_xtnd effective diffusion coefficient from cg_drag_mod ! [ m^2/s ] ! source_level k index of gravity wave source 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 ! frequency [ /s^2 ] ! ie, je ending subdomain indices of data in the current ! physics window being integrated ! imax, jmax, kmax ! physics window dimensions ! i, j, k, nn do loop indices ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! define processor extents and loop limits. !--------------------------------------------------------------------- imax = size(uuu,1) jmax = size(uuu,2) kmax = size(uuu,3) ie = is + imax - 1 je = js + jmax - 1 !--------------------------------------------------------------------- ! 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 !----------------------------------------------------------------------- ! 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) dtdz(i,j,1) = (temp (i,j,1) - temp (i,j,2))/ & (zfull(i,j,1) - zfull(i,j,2)) do k=2,iz0 dtdz(i,j,k) = (temp (i,j,k-1) - temp (i,j,k+1))/ & (zfull(i,j,k-1) - zfull(i,j,k+1)) end do !-------------------------------------------------------------------- ! calculate air density. !-------------------------------------------------------------------- do k=1,iz0+1 zden(i,j,k ) = pfull(i,j,k)/(temp(i,j,k)*RDGAS) end do !---------------------------------------------------------------------- ! calculate buoyancy frequency. restrict the squared buoyancy ! frequency to be no smaller than bflim. !---------------------------------------------------------------------- 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 zbf(i,j,k) = sqrt(zbf(i,j,k)) endif end do !---------------------------------------------------------------------- ! if zbf is to be saved for netcdf output, the remaining vertical ! levels must be initialized. !---------------------------------------------------------------------- if (id_bf_cgwd > 0) then zbf(i,j,iz0+1:) = 0.0 endif !---------------------------------------------------------------------- ! define an array of heights at model levels and an array containing ! the zonal wind component. !---------------------------------------------------------------------- do k=1,iz0+1 zzchm(i,j,k) = zfull(i,j,k) end do do k=1,iz0 zu(i,j,k) = uuu(i,j,k) zv(i,j,k) = vvv(i,j,k) end do !---------------------------------------------------------------------- ! 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 ! frequency - constant across model top. !---------------------------------------------------------------------- zzchm(i,j,0) = zzchm(i,j,1) + zzchm(i,j,1) - zzchm(i,j,2) zu(i,j,0) = 2.*zu(i,j,1) - zu(i,j,2) zv(i,j,0) = 2.*zv(i,j,1) - zv(i,j,2) zden(i,j,0) = zden(i,j,1)*zden(i,j,1)/zden(i,j,2) 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 ! 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. !--------------------------------------------------------------------- call gwfc (is, ie, js, je, source_level, source_amp, & 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, & 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. !------------------------------------------------------------------- gwd_u(is:ie,js:je,:) = gwfcng_x(:,:,:) 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 !-------------------------------------------------------------------- ! 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) endif if (id_kedy_cgwd > 0) then used = send_data (id_kedy_cgwd, ked_gwfc_y, Time, is, js, 1) 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 ) endif if (id_gwfx_cgwd > 0) then 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) endif !-------------------------------------------------------------------- ! 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) gwfcng_x(:,:,:) = gwd_u(is:ie,js:je,:) gwfcng_y(:,:,:) = gwd_v(is:ie,js:je,:) endif ! (cgdrag_alarm <= 0) !-------------------------------------------------------------------- end subroutine cg_drag_calc !################################################################### subroutine cg_drag_end !-------------------------------------------------------------------- ! cg_drag_end is the destructor for cg_drag_mod. !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables !For version 3 and after, use NetCDF restarts. if (mpp_pe() == mpp_root_pe() ) & call error_mesg ('cg_drag_mod', 'write_restart_nc: & &Writing netCDF formatted restart file as & &requested. ', NOTE) call cg_drag_restart #ifdef COL_DIAG if (column_diagnostics_desired) then call close_column_diagnostics_units (diag_units) endif #endif !--------------------------------------------------------------------- ! mark the module as uninitialized. !--------------------------------------------------------------------- module_is_initialized = .false. !--------------------------------------------------------------------- end subroutine cg_drag_end !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! PRIVATE SUBROUTINES ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine read_nc_restart_file !----------------------------------------------------------------------- ! subroutine read_restart_nc reads a netcdf restart file to obtain ! the variables needed upon experiment restart. !----------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables: character(len=64) :: fname='INPUT/cg_drag.res.nc' character(len=8) :: chvers !--------------------------------------------------------------------- ! local variables: ! ! fname restart file name ! !---------------------------------------------------------------------- !-------------------------------------------------------------------- ! output a message indicating entrance into this routine. !-------------------------------------------------------------------- if (mpp_pe() == mpp_root_pe() ) then call error_mesg ('cg_drag_mod', 'read_restart_nc:& &Reading netCDF formatted restart file:'//trim(fname), NOTE) endif !------------------------------------------------------------------- ! read the values of gwd_u and gwd_v !------------------------------------------------------------------- if (size(restart_versions(:)) .le. 2 ) then call error_mesg ('cg_drag_mod', 'read_restart_nc: restart file format is netcdf, ' // & 'restart_versions is not netcdf file version', FATAL) endif call restore_state(Cg_restart) if(in_different_file) call restore_state(Til_restart) if (.not. any(vers == restart_versions) ) then write (chvers, '(i4)') vers call error_mesg ('cg_drag_init', & 'restart version '//chvers//' cannot be read & &by this module version', FATAL) endif 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 (cg_drag_freq /= old_time_step) then cgdrag_alarm = cgdrag_alarm - old_time_step + cg_drag_freq if (mpp_pe() == mpp_root_pe() ) then call error_mesg ('cg_drag_mod', & 'cgdrag time step has changed, & &next cgdrag time also changed', NOTE) endif old_time_step = cg_drag_freq endif !-------------------------------------------------------------------- ! if cg_drag_offset is specified and is smaller than the time remain- ! 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. !-------------------------------------------------------------------- if (cg_drag_offset /= 0) then if (cgdrag_alarm > cg_drag_offset) then cgdrag_alarm = cg_drag_offset endif endif !--------------------------------------------------------------------- end subroutine read_nc_restart_file !#################################################################### ! register restart field to be read and written through save_restart and restore_state. subroutine cg_drag_register_restart character(len=64) :: fname = 'cg_drag.res.nc' ! name of restart file character(len=64) :: fname2 integer :: id_restart call get_mosaic_tile_file(fname, fname2, .false. ) allocate(Cg_restart) if(trim(fname2) == trim(fname)) then Til_restart => Cg_restart in_different_file = .false. else in_different_file = .true. allocate(Til_restart) endif id_restart = register_restart_field(Cg_restart, fname, 'restart_version', vers, no_domain = .true. ) id_restart = register_restart_field(Cg_restart, fname, 'cgdrag_alarm', cgdrag_alarm, no_domain = .true. ) id_restart = register_restart_field(Cg_restart, fname, 'cg_drag_freq', old_time_step, no_domain = .true. ) id_restart = register_restart_field(Til_restart, fname, 'gwd_u', gwd_u) id_restart = register_restart_field(Til_restart, fname, 'gwd_v', gwd_v) return end subroutine cg_drag_register_restart !#################################################################### ! ! ! ! write out restart file. ! 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. ! ! subroutine cg_drag_restart(timestamp) character(len=*), intent(in), optional :: timestamp call save_restart(Cg_restart, timestamp) if(in_different_file) call save_restart(Til_restart, timestamp) end subroutine cg_drag_restart ! NAME=cg_drag_restart" !#################################################################### subroutine gwfc (is, ie, js, je, source_level, source_amp, 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. ! Alexander and L. Bruhwiler. !------------------------------------------------------------------- !------------------------------------------------------------------- integer, intent(in) :: is, ie, js, je integer, dimension(:,:), intent(in) :: source_level real, dimension(:,:), intent(in) :: source_amp real, dimension(:,:,0:), intent(in) :: rho, u, bf, z real, dimension(:,:,0:), intent(out) :: gwf real, dimension(:,:,0:), intent(out) :: ked !------------------------------------------------------------------- ! intent(in) variables: ! ! is, ie, js, je starting/ending subdomain i,j indices of data ! in the physics_window being integrated ! source_level k index of model level serving as gravity wave ! source ! source_amp amplitude of gravity wave source ! ! rho atmospheric density [ kg/m^3 ] ! u zonal wind component [ m/s ] ! bf buoyancy frequency [ /s ] ! z height of model levels [ m ] ! ! intent(out) variables: ! ! gwf gravity wave forcing in u equation [ m/s^2 ] ! ! intent(out), optional variables: ! ! ked eddy diffusion coefficient from gravity wave ! forcing [ m^2/s ] ! !------------------------------------------------------------------ !------------------------------------------------------------------ ! local variables real, dimension (0:size(u,3)-1 ) :: & wv_frcng, diff_coeff, c0mu, dz, & fac, omc integer, dimension (nc) :: msk real , dimension (nc) :: c0mu0, B0 real :: fm, fe, Hb, alp2, Foc, c, test, rbh,& eps, Bsum integer :: iz0 integer :: i, j, k, ink, n real :: ampl !------------------------------------------------------------------ ! 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 ! [ m/s ] ! dz delta z between model levels [ m ] ! fac factor used in determining if wave is breaking ! [ s/m ] ! 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 ! 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 ! 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 ] ! Hb density scale height [ m ] ! alp2 scale height factor: 1/(2*Hb)**2 [ /m^2 ] ! Foc wave breaking threshold [ s/m ] ! c wave phase speed used in defining wave momentum flux ! amplitude [ m/s ] ! test condition defining internal reflection [ /s ] ! rbh atmospheric density at half-level (geometric mean) ! [ kg/m^3 ] ! eps intermittency factor ! 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 phase speed loop index ! !-------------------------------------------------------------------- !------------------------------------------------------------------- ! 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. !------------------------------------------------------------------- gwf = 0.0 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) !-------------------------------------------------------------------- ! 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) !--------------------------------------------------------------------- ! when the wave phase speed is same as wind speed, there is no ! momentum flux. !--------------------------------------------------------------------- if (c0mu0(n) == 0.0) then B0(n) = 0.0 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) 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)) endif Bsum = Bsum + abs(B0(n)) endif end do !--------------------------------------------------------------------- ! define the intermittency factor eps. the factor of 1.5 is currently ! unexplained. !--------------------------------------------------------------------- if (Bsum == 0.0) then call error_mesg ('cg_drag_mod', & ' zero flux input at source level', FATAL) endif eps = (ampl*1.5/nk)/Bsum !-------------------------------------------------------------------- ! loop over the nk different wavelengths in the spectrum. !-------------------------------------------------------------------- do ink=1,nk ! wavelength loop !---------------------------------------------------------------------- ! define variables needed at levels above the source level. !--------------------------------------------------------------------- do k=0,iz0 fac(k) = 0.5*(rho(i,j,k)/rho(i,j,iz0))*kwv(ink)/bf(i,j,k) end do 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) omc(k) = sqrt((bf(i,j,k)*bf(i,j,k)*k2(ink))/ & (k2(ink) + alp2)) end do !--------------------------------------------------------------------- ! 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 ! from all waves breaking at a given level. !---------------------------------------------------------------------- do k=iz0, 0, -1 fm = 0. fe = 0. do n=1,nc ! phase speed loop !---------------------------------------------------------------------- ! check only those waves which are still propagating, i.e., msk = 1. !---------------------------------------------------------------------- if (msk(n) == 1) then c0mu(k) = c0(n) - u(i,j,k) !---------------------------------------------------------------------- ! if phase speed matches the wind speed, remove c0(n) from the ! set of propagating waves. !---------------------------------------------------------------------- if (c0mu(k) == 0.) then msk(n) = 0 else !--------------------------------------------------------------------- ! define the criterion which determines if wave is reflected at this ! level (test). !--------------------------------------------------------------------- test = abs(c0mu(k))*kwv(ink) - omc(k) if (test >= 0.0) then !--------------------------------------------------------------------- ! wave has undergone total internal reflection. remove it from the ! propagating set. !--------------------------------------------------------------------- msk(n) = 0 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 ! (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. ! set flag to remove phase speed c0(n) from the set of active waves ! moving upwards to the next level. !--------------------------------------------------------------------- Foc = B0(n)/(c0mu(k) )**3 - fac(k) if ((Foc >= 0.0) .or. & (c0mu0(n)*c0mu(k) <= 0.0)) then msk(n) = 0 if (k < iz0) then fm = fm + B0(n) fe = fe + c0mu(k)*B0(n) endif endif endif ! (test >= 0.0) endif ! (c0mu == 0.0) endif ! (msk == 1) end do ! phase speed loop !---------------------------------------------------------------------- ! 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) wv_frcng(k+1) = 0.5*(wv_frcng(k+1) + wv_frcng(k)) diff_coeff(k) = (rho(i,j,iz0)/rbh)*fe*eps/(dz(k)* & bf(i,j,k)*bf(i,j,k)) diff_coeff(k+1) = 0.5*(diff_coeff(k+1) + & diff_coeff(k)) else wv_frcng(iz0) = 0.0 diff_coeff(iz0) = 0.0 endif end do ! (k loop) !--------------------------------------------------------------------- ! increment the total forcing at each point with that obtained from ! the set of waves with the current wavenumber. !--------------------------------------------------------------------- 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 ! wavelength loop end do ! i loop end do ! j loop !-------------------------------------------------------------------- end subroutine gwfc !#################################################################### end module cg_drag_mod