Merge pull request #9 from mjucker/local_heating

Local heating

.gitignore

@@ -29,7 +29,6 @@ git_message
 *~
 it_jobscript*
 simple_jobscript*
-local_*
 bin/mppnccombine.nyu
 
 .DS_Store

docs/Parameters.md

@@ -227,6 +227,7 @@ The namelist `physics_driver_nml` steers which physics components are used. It r
  do_grey_radiation | .false. | rather do grey radiation?
  do_rrtm_radiation | .true. | or RRTM radiation?
  do_damping | .true. | do any of the damping schemes?
+ do_local_heating | .false. | add artificial local heating? If so, see `local_heating_nml` namelist
  diff_min | 1.e-3    | minimum value of a diffusion coefficient beneath which the coefficient is reset to zero
  diffusion_smooth | .true. | diffusion coefficients should be smoothed in time?
  do_netcdf_restart | .true. | make restart files netCDF format?
@@ -437,6 +438,20 @@ Variable | Default Value | Meaning
         equinox_day | 0.25              | fraction of the year defining March equinox.
 
 
+#### Local heating
+      
+
+If `do_local_heating = .true.` in `physics_driver_nml`, the namelist `local_heating_nml` can be used to set the form and position of the desired local heating.
+
+Variable | Default Value | Meaning
+ :--- | :---: | :---
+      hamp | 0 | amplitude of Gaussian heating in [K/s], maximum 10 entries
+      latcenter| 90 | horizontal center of the Gaussian in [degrees latitude], maximum 10 entries
+      latwidth | 15 | horizontal width of the Gaussian in the [degrees latitude], maximum 10 entries
+      pcenter | 1 | vertical center of the Gaussian in the vertical [hPa], maximum 10 entries
+      pwidth  | 2  | vertical width of the Gaussian in orders of magnitude [log10(hPa)], maximum 10 entries
+
+
 ### Boundary conditions
 
 The upper boundary conditions are defined by how to deal with gravity waves. The simplest choice is simple Rayleigh friction (damping towards zero momentum), non-orographic gravity wave parameterization, orographic gravity wave drag, or a constant "gravity wave" drag.

exp/path_names

@@ -1,5 +1,6 @@
 atmos_coupled/atmos_model.f90
 atmos_param/qflux/qflux.f90
+atmos_param/local_heating/local_heating.f90 
 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.f90
 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parrrtm.f90
 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/rrlw_cld.f90

src/atmos_param/local_heating/local_heating.f90

@@ -0,0 +1,167 @@
+module local_heating_mod
+
+
+  use           fms_mod, only: error_mesg, FATAL, file_exist,       &
+                               open_namelist_file, set_domain,      &
+  			                       read_data, check_nml_error,          &
+                               mpp_pe, mpp_root_pe, close_file,     &
+                               write_version_number, stdlog,        &
+                               uppercase,&  !pjk
+                               mpp_clock_id,mpp_clock_begin,mpp_clock_end,CLOCK_COMPONENT!,mpp_chksum
+
+  use diag_manager_mod, only: register_diag_field, send_data
+  use time_manager_mod, only: time_type
+
+  use  time_manager_mod, only: time_type, get_time, length_of_year
+
+  use  diag_manager_mod, only: register_diag_field, send_data
+
+  use  field_manager_mod, only: MODEL_ATMOS, parse
+
+  use constants_mod, only   : RADIAN,PI
+
+  use rrtm_astro,only : equinox_day
+
+  implicit none
+
+  !-----------------------------------------------------------------------
+  !---------- interfaces ------------
+  public :: local_heating,local_heating_init
+
+  !---------------------------------------------------------------------------------------------------------------
+  !
+  !-------------------- diagnostics fields -------------------------------
+  integer :: id_tdt_lheat
+  character(len=14) :: mod_name = 'local_heating'
+  real :: missing_value = -999.
+  !-----------------------------------------------------------------------
+  !-------------------- namelist -----------------------------------------
+  !-----------------------------------------------------------------------
+  integer,parameter :: ngauss = 10
+  real,dimension(ngauss)   :: hamp      = 0.        ! heating amplitude [K/d]
+  real,dimension(ngauss)   :: latwidth  = 15.       ! latitudinal width of Gaussian heating [deg]
+  real,dimension(ngauss)   :: latcenter = 90.       ! latitudinal center of Gaussian heating [deg]
+  real,dimension(ngauss)   :: pwidth    = 2.        ! height of Gaussian heating in log-pressure [log10(hPa)]
+  real,dimension(ngauss)   :: pcenter   = 1.        ! center of Gaussian heating in pressure [hPa]
+  integer,dimension(ngauss):: hk        = 0         ! temporal wave number for cosine
+                                                    !  0 -> no cosine
+                                                    !  1 -> heating in one season only
+                                                    !  2 -> heating in two seasons, etc
+  real,dimension(ngauss)   :: hphase    = 0.        ! phasing of cosine in annual cycle in [pi]
+                                                    !  relative to March equinox
+                                                    !  example: for hk = 1, that means
+                                                    !  0   -> heating in JFMAMJ
+                                                    !  0.5 -> heating in AMJJAS
+                                                    !  1   -> heating in JASOND
+                                                    !  1.5 -> heating in ONDJFM
+
+  namelist /local_heating_nml/ hamp,latwidth,latcenter,pwidth,pcenter,hk,hphase
+
+
+  ! local variables
+  real,dimension(ngauss) :: logpc
+
+contains
+
+  subroutine local_heating_init(axes, Time)
+    implicit none
+    integer, intent(in), dimension(4) :: axes
+    type(time_type), intent(in)       :: Time
+    !-----------------------------------------------------------------------
+    integer :: unit, io, ierr, n
+
+    !     ----- read namelist -----
+
+    if (file_exist('input.nml')) then
+       unit = open_namelist_file ( )
+       ierr=1
+       do while (ierr /= 0)
+          read  (unit, nml=local_heating_nml, iostat=io, end=10)
+          ierr = check_nml_error (io, 'local_heating_nml')
+       enddo
+10     call close_file (unit)
+    endif
+
+   ! ---- convert input units to code units  -----
+    do n = 1,ngauss
+       pcenter(n)   = pcenter(n)*100      ! convert hPa to Pa
+       logpc(n)     = log10(pcenter(n))   ! work in log10(p)
+       hamp(n)      = hamp(n)/86400.      ! convert K/d to K/s
+       latcenter(n) = latcenter(n)/RADIAN ! convert degrees to radians
+       latwidth(n)  = latwidth(n)/RADIAN  ! convert degrees to radians
+       hphase(n)    = hphase(n)*PI        ! convert to radians
+    enddo
+  !----
+  !------------ initialize diagnostic fields ---------------
+    id_tdt_lheat = &
+         register_diag_field ( mod_name, 'tdt_lheat', axes(1:3), Time, &
+         'Temperature tendency due to local heating', &
+         'K/s', missing_value=missing_value               )
+
+  end subroutine local_heating_init
+
+  !-----------------------------------------------------------------------
+  !-------------------- computing localized heating ----------------------
+  !-----------------------------------------------------------------------
+
+  subroutine local_heating(is,js,Time,lon,lat,p_full,tdt_tot)
+    implicit none
+    integer, intent(in)                  :: is, js
+    type(time_type),intent(in)           :: Time
+    real, dimension(:,:)  ,intent(in)    :: lon,lat
+    real, dimension(:,:,:),intent(in)    :: p_full
+    real, dimension(:,:,:),intent(inout) :: tdt_tot
+    ! local variables
+    integer :: i,j,k,n
+    real, dimension(size(lon,1),size(lon,2)) :: lat_factor
+    real, dimension(size(tdt_tot,1),size(tdt_tot,2),size(tdt_tot,3)) :: tdt
+    real    :: logp,p_factor,t_factor,cosdays
+    logical :: used
+    integer :: seconds,days,daysperyear
+
+    ! get temporal dependencies
+    call get_time(length_of_year(),seconds,daysperyear)
+    call get_time(Time,seconds,days)
+    !  day of the year relative to equinox_day
+    days = days - int(equinox_day*daysperyear)
+    cosdays = modulo(days,daysperyear)*1.0/daysperyear*2*PI
+    ! days is now the day of the year relative to March equinox
+
+    tdt = 0.
+    do n = 1,ngauss
+       if ( hamp(n) .ne. 0. ) then
+          ! temporal component
+          t_factor = max(0.0,cos(hk(n)*cosdays-hphase(n)))
+          ! meridional component
+          do j = 1,size(lon,2)
+             do i = 1,size(lon,1)
+                lat_factor(i,j) = exp( -(lat(i,j)-latcenter(n))**2/(2*(latwidth(n))**2) )
+             enddo
+          enddo
+
+          do k=1,size(p_full,3)
+             do j = 1,size(lon,2)
+                do i = 1,size(lon,1)
+                   logp = log10(p_full(i,j,k))
+                   ! vertical component
+                   p_factor = exp(-(logp-logpc(n))**2/(2*(pwidth(n))**2))
+                   ! everything together
+                   tdt(i,j,k) = tdt(i,j,k) + hamp(n)*t_factor*lat_factor(i,j)*p_factor
+                enddo
+             enddo
+          enddo
+       endif
+    enddo
+
+    tdt_tot = tdt_tot + tdt
+
+   !------- diagnostics ------------
+    if ( id_tdt_lheat > 0 ) then
+       used = send_data ( id_tdt_lheat, tdt, Time, is, js, 1 )
+    endif
+
+   end subroutine local_heating
+  !-----------------------------------------------------------------------
+  !-----------------------------------------------------------------------
+
+end module local_heating_mod

src/atmos_param/physics_driver/physics_driver.f90

@@ -126,6 +126,8 @@ module physics_driver_mod
 
 use grey_radiation_mod, only: grey_radiation_init, grey_radiation, grey_radiation_end
 
+use local_heating_mod, only:  local_heating_init,local_heating
+
 !mj: RRTM radiative scheme
 use rrtmg_lw_init
 use rrtmg_lw_rad
@@ -193,6 +195,8 @@ module physics_driver_mod
 
 logical :: do_damping = .true.
 
+logical :: do_local_heating = .false.
+
 real    :: diff_min = 1.e-3    ! minimum value of a diffusion 
                                ! coefficient beneath which the
                                ! coefficient is reset to zero
@@ -230,7 +234,7 @@ module physics_driver_mod
                                 do_moist_processes, tau_diff,      &
                                 diff_min, diffusion_smooth, &
                                 do_grey_radiation, do_rrtm_radiation, &
-                                do_damping
+                                do_damping, do_local_heating
 
 !---------------------------------------------------------------------
 !------- public data ------
@@ -573,6 +577,8 @@ subroutine physics_driver_init (Time, lonb, latb, axes, pref, &
          call rrtmg_sw_ini(cp_air)
          call rrtm_radiation_init(axes,Time,id*jd,kd,lonb,latb)
       endif
+
+      if(do_local_heating) call local_heating_init(axes, Time)
 
 !-----------------------------------------------------------------------
 !    initialize atmos_tracer_driver_mod.
@@ -1282,6 +1288,12 @@ subroutine physics_driver_down (is, ie, js, je,                       &
          call interp_temp(z_full,z_half,t_surf_rad,t)
          call run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt,coszen,flux_sw,flux_lw)
       endif
+!----------------------------------------------------------------------
+!    artificial local heating if required
+!----------------------------------------------------------------------
+      if(do_local_heating) then
+        call local_heating(is,js,Time,lon,lat,p_full,tdt)
+      endif
 
 !----------------------------------------------------------------------
 !    call damping_driver to calculate the various model dampings that