mixed_layer_mod.F90

module mixed_layer_mod

!
! Implementation of mixed layer boundary condition
!

!use                  fms_mod, only: set_domain, write_version_number, &
!                                    mpp_pe, mpp_root_pe, error_mesg, FATAL, WARNING

!use                  fms_mod, only: stdlog, check_nml_error, close_file,&
!                                    open_namelist_file, stdout, file_exist, &
!                                    read_data, write_data, open_file

use           gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit
use            constants_mod, only: HLV, PI, RHO_CP, CP_AIR

!use         diag_manager_mod, only: register_diag_field, send_data

!use       time_manager_mod,   only: time_type

!use           transforms_mod, only: get_deg_lat

!use            vert_diff_mod, only: surf_diff_type

implicit none
private
!===================================================================================================

character(len=128) :: version = '$Id: mixed_layer_mod.F90,v 1.3 2017/08/11 20:48:51 jmc Exp $'
character(len=128) :: tagname = '$Name:  $'
character(len=128), parameter :: mod_name='mixed_layer'

!===================================================================================================

public :: mixed_layer_init, mixed_layer, mixed_layer_end

!===================================================================================================

logical :: evaporation = .true.
real    :: qflux_amp = 0.0
real    :: qflux_width = 16.0  ! width of qflux region in degrees
real    :: depth = 40.0

namelist/mixed_layer_nml/ evaporation, qflux_amp, depth, qflux_width

!===================================================================================================

logical :: module_is_initialized =.false.
logical :: used

integer :: iter
integer, dimension(4) :: axes
!integer ::                                                                    &
!     id_t_surf,            &   ! surface temperature
!     id_flux_lhe,          &   ! latent heat flux at surface
!     id_flux_oceanq,       &   ! oceanic Q flux
!     id_flux_t                 ! sensible heat flux at surface

!real, allocatable, dimension(:,:)   ::                                        &
!     ocean_qflux,           &   ! Q-flux
!     rad_lat_2d                 ! latitude in radians

!real, allocatable, dimension(:)   :: deg_lat

!real, allocatable, dimension(:,:)   ::                                        &
!     gamma_t,               &   ! Used to calculate the implicit
!     gamma_q,               &   ! correction to the diffusion in
!     fn_t,                  &   ! the lowest layer
!     fn_q,                  &   !
!     en_t,                  &   !
!     en_q,                  &   !
!     alpha_t,               &   !
!     alpha_q,               &   !
!     alpha_lw,              &   !
!     beta_t,                &   !
!     beta_q,                &   !
!     beta_lw,               &   !
!     t_surf_dependence,     &   !
!     corrected_flux,        &   !
!     eff_heat_capacity,     &   ! Effective heat capacity
!     delta_t_surf               ! Increment in surface temperature

real inv_cp_air

!===================================================================================================
contains
!===================================================================================================

!subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, axes,  &
!                            rad_lat_2d, ocean_qflux,                   &
subroutine mixed_layer_init( is, ie, js, je, num_levels,         axes,  &
                             Time, cst_mxlDepth, myThid )

integer, intent(in) :: is, ie, js, je, num_levels
!real, intent(inout), dimension(:,:) :: t_surf
integer, intent(in), dimension(4) :: axes
!type(time_type), intent(in)       :: Time
real, intent(in)                   :: Time
real, intent(out)                  :: cst_mxlDepth
integer, intent(in)                :: myThid
!real, intent(in), dimension(:,:)  :: rad_lat_2d
!real, intent(out), dimension(:,:) :: ocean_qflux

!integer :: j
!real    :: rad_qwidth
!integer:: ierr, io, unit
integer         :: iUnit
CHARACTER*(gcm_LEN_MBUF) :: msgBuf

if(module_is_initialized) return

! read namelist and copy to logfile
!    _BARRIER
!    _BEGIN_MASTER(myThid)
     CALL BARRIER(myThid)
     IF ( myThid.EQ.1 ) THEN

     WRITE(msgBuf,'(A)') 'MIXED_LAYER_INIT: opening data.atm_gray'
     CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
     CALL OPEN_COPY_DATA_FILE(                                      &
                           'data.atm_gray', 'MIXED_LAYER_INIT',     &
                           iUnit,                                   &
                           myThid )
!    Read parameters from open data file
     READ(UNIT=iUnit,NML=mixed_layer_nml)
     WRITE(msgBuf,'(A)')                                            &
          'MIXED_LAYER_INIT: finished reading data.atm_gray'
     CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
!    Close the open data file
#ifdef SINGLE_DISK_IO
     CLOSE(iUnit)
#else
     CLOSE(iUnit,STATUS='DELETE')
#endif /* SINGLE_DISK_IO */

!call write_version_number(version, tagname)

!unit = open_namelist_file ()
!ierr=1
!do while (ierr /= 0)
!  read  (unit, nml=mixed_layer_nml, iostat=io, end=10)
!  ierr = check_nml_error (io, 'mixed_layer_nml')
!enddo
!10 call close_file (unit)

!if ( mpp_pe() == mpp_root_pe() )   write (stdlog(), nml=mixed_layer_nml)

!allocate(rad_lat_2d              (is:ie, js:je))
!allocate(ocean_qflux             (is:ie, js:je))

!allocate(deg_lat                 (js:je))

!allocate(gamma_t                 (is:ie, js:je))
!allocate(gamma_q                 (is:ie, js:je))
!allocate(en_t                    (is:ie, js:je))
!allocate(en_q                    (is:ie, js:je))
!allocate(fn_t                    (is:ie, js:je))
!allocate(fn_q                    (is:ie, js:je))
!allocate(alpha_t                 (is:ie, js:je))
!allocate(alpha_q                 (is:ie, js:je))
!allocate(alpha_lw                (is:ie, js:je))
!allocate(beta_t                  (is:ie, js:je))
!allocate(beta_q                  (is:ie, js:je))
!allocate(beta_lw                 (is:ie, js:je))
!allocate(delta_t_surf            (is:ie, js:je))
!allocate(eff_heat_capacity       (is:ie, js:je))
!allocate(corrected_flux          (is:ie, js:je))
!allocate(t_surf_dependence       (is:ie, js:je))

inv_cp_air = 1.0 / CP_AIR

module_is_initialized = .true.

     ENDIF
     CALL BARRIER(myThid)
     cst_mxlDepth = depth

!see if restart file exists for the surface temperature
!
!if (file_exist('INPUT/mixed_layer.res')) then
!         unit = open_file (file='INPUT/mixed_layer.res', &
!                           form='native', action='read')
!         call read_data (unit, t_surf)
!         call close_file (unit)
!else if (file_exist('INPUT/swamp.res')) then
!         unit = open_file (file='INPUT/swamp.res', &
!                           form='native', action='read')
!         call read_data (unit, t_surf)
!         call close_file (unit)
!  call error_mesg('mixed_layer','mixed_layer restart file not found, using swamp restart file', WARNING)
!else
!  call error_mesg('mixed_layer','mixed_layer restart file not found', WARNING)
!endif

!id_t_surf = register_diag_field(mod_name, 't_surf',        &
!                                axes(1:2), Time, 'surface temperature','K')
!id_flux_t = register_diag_field(mod_name, 'flux_t',        &
!                                axes(1:2), Time, 'sensible heat flux up at surface','watts/m2')
!id_flux_lhe = register_diag_field(mod_name, 'flux_lhe',        &
!                                 axes(1:2), Time, 'latent heat flux up at surface','watts/m2')
!id_flux_oceanq = register_diag_field(mod_name, 'flux_oceanq',        &
!                                 axes(1:2), Time, 'oceanic Q-flux','watts/m2')

! latitude will be needed for oceanic q flux
!call get_deg_lat(deg_lat)
!do j=js,je
!  rad_lat_2d(:,j) = deg_lat(j)*PI/180.
!enddo

! calculate ocean Q flux
!rad_qwidth = qflux_width*PI/180.
!ocean_qflux = qflux_amp*(1-2.*rad_lat_2d**2/rad_qwidth**2) * &
!        exp(- ((rad_lat_2d)**2/(rad_qwidth)**2))

return
end subroutine mixed_layer_init

!===================================================================================================

subroutine mixed_layer (                                &
     Time,                                              &
     t_surf,                                            &
     flux_t,                                            &
     flux_q,                                            &
     flux_r,                                            &
     dt,                                                &
     net_surf_sw_down,                                  &
     surf_lw_down,                                      &
!    Tri_surf,                                          &
                tri_surf_dtmass,                               &
                tri_surf_dflux_t, tri_surf_dflux_q,            &
                tri_surf_delta_t, tri_surf_delta_q,            &
     dhdt_surf,                                         &
     dedt_surf,                                         &
     dedq_surf,                                         &
     drdt_surf,                                         &
     dhdt_atm,                                          &
     dedq_atm,                                          &
     ocean_qflux, mixLayDepth,                          &
     delta_t_surf, &               ! Increment in surface temperature
     myThid )

! ---- arguments -----------------------------------------------------------
!type(time_type), intent(in)       :: Time
real, intent(in)                  :: Time
real, intent(in),  dimension(:,:) :: &
     net_surf_sw_down, surf_lw_down
real, intent(inout), dimension(:,:) :: &
     flux_t,    flux_q,     flux_r
real, intent(inout), dimension(:,:) :: t_surf
real, intent(in), dimension(:,:) :: &
   dhdt_surf, dedt_surf, dedq_surf, &
   drdt_surf, dhdt_atm, dedq_atm
real, intent(in) :: dt
!type(surf_diff_type), intent(inout) :: Tri_surf
real, intent(inout), dimension(:,:)  :: tri_surf_dtmass
real, intent(inout), dimension(:,:)  :: tri_surf_dflux_t, tri_surf_dflux_q
real, intent(inout), dimension(:,:)  :: tri_surf_delta_t, tri_surf_delta_q
real, intent(in),    dimension(:,:)  :: ocean_qflux
real, intent(in),    dimension(:,:)  :: mixLayDepth
real, intent(out),   dimension(:,:)  :: delta_t_surf
integer, intent(in) :: myThid

! ---- local variables --------------------------------------------------
real, allocatable, dimension(:,:)   ::                                        &
     gamma_t,               &   ! Used to calculate the implicit
     gamma_q,               &   ! correction to the diffusion in
     fn_t,                  &   ! the lowest layer
     fn_q,                  &   !
     en_t,                  &   !
     en_q,                  &   !
     alpha_t,               &   !
     alpha_q,               &   !
     alpha_lw,              &   !
     beta_t,                &   !
     beta_q,                &   !
     beta_lw,               &   !
     t_surf_dependence,     &   !
     corrected_flux,        &   !
     eff_heat_capacity          ! Effective heat capacity
!    eff_heat_capacity,     &   ! Effective heat capacity
!    delta_t_surf               ! Increment in surface temperature
integer :: im, jm
!----------------------------------------------------------------

im = size(t_surf,1)
jm = size(t_surf,2)
allocate(gamma_t            (im,jm) )
allocate(gamma_q            (im,jm) )
allocate(en_t               (im,jm) )
allocate(en_q               (im,jm) )
allocate(fn_t               (im,jm) )
allocate(fn_q               (im,jm) )
allocate(alpha_t            (im,jm) )
allocate(alpha_q            (im,jm) )
allocate(alpha_lw           (im,jm) )
allocate(beta_t             (im,jm) )
allocate(beta_q             (im,jm) )
allocate(beta_lw            (im,jm) )
allocate(eff_heat_capacity  (im,jm) )
allocate(corrected_flux     (im,jm) )
allocate(t_surf_dependence  (im,jm) )
!allocate(delta_t_surf       (im,jm) )

if(.not.module_is_initialized) then
! call error_mesg('mixed_layer','mixed_layer module is not initialized',FATAL)
  call PRINT_ERROR('mixed_layer: mixed_layer module is not initialized',myThid)
  STOP 'ABNORMAL END: S/R MIXED_LAYER'
endif

! Need to calculate the implicit changes to the lowest level delta_q and delta_t
! - see the discussion in vert_diff.tech.ps

! Care is needed to differentiate between the sensible heat flux and the
! diffusive flux of temperature

gamma_t = 1.0 / (1.0 - Tri_surf_dtmass * (Tri_surf_dflux_t + dhdt_atm * inv_cp_air))
gamma_q = 1.0 / (1.0 - Tri_surf_dtmass * (Tri_surf_dflux_q + dedq_atm))

fn_t = gamma_t * (Tri_surf_delta_t + Tri_surf_dtmass * flux_t * inv_cp_air)
fn_q = gamma_q * (Tri_surf_delta_q + Tri_surf_dtmass * flux_q)

en_t = gamma_t * Tri_surf_dtmass * dhdt_surf * inv_cp_air
en_q = gamma_q * Tri_surf_dtmass * dedt_surf

!
! Note flux_sw doesn't depend on surface or lowest layer values
! Note drdt_atm is not used - should be fixed
!
alpha_t = flux_t * inv_cp_air + dhdt_atm * inv_cp_air * fn_t
alpha_q = flux_q + dedq_atm * fn_q
alpha_lw = flux_r

beta_t = dhdt_surf * inv_cp_air + dhdt_atm * inv_cp_air * en_t
beta_q = dedt_surf + dedq_atm * en_q
beta_lw = drdt_surf

!
! Implement mixed layer surface boundary condition
!
corrected_flux = - net_surf_sw_down - surf_lw_down + alpha_t * CP_AIR + alpha_lw + ocean_qflux
t_surf_dependence = beta_t * CP_AIR + beta_lw

if (evaporation) then
  corrected_flux = corrected_flux + alpha_q * HLV
  t_surf_dependence = t_surf_dependence + beta_q * HLV
endif

!
! Now update the mixed layer surface temperature using an implicit step
!
!eff_heat_capacity = depth * RHO_CP + t_surf_dependence * dt
eff_heat_capacity = mixLayDepth * RHO_CP + t_surf_dependence * dt

if (any(eff_heat_capacity .eq. 0.0))  then
  write(*,*) 'mixed_layer: error', eff_heat_capacity
! call error_mesg('mixed_layer', 'Avoiding division by zero',fatal)
  call PRINT_ERROR('mixed_layer: Avoiding division by zero',myThid)
  STOP 'ABNORMAL END: S/R MIXED_LAYER'
end if

delta_t_surf = - corrected_flux  * dt / eff_heat_capacity

t_surf = t_surf + delta_t_surf

! Finally calculate the increments for the lowest atmospheric layer
!
Tri_surf_delta_t = fn_t + en_t * delta_t_surf
Tri_surf_delta_q = fn_q + en_q * delta_t_surf

! Note:
! When using an implicit step there is not a clearly defined flux for a given timestep

! Final flux values (match mixed layer heat content variations)
  flux_t = ( alpha_t  + delta_t_surf*beta_t )*CP_AIR
  flux_r =   alpha_lw + delta_t_surf*beta_lw
  flux_q =   alpha_q  + delta_t_surf*beta_q

!if(id_t_surf > 0) used = send_data(id_t_surf, t_surf, Time)
!if(id_flux_t > 0) used = send_data(id_flux_t, flux_t, Time)
!if(id_flux_lhe > 0) used = send_data(id_flux_lhe, HLV * flux_q, Time)
!if(id_flux_oceanq > 0)   used = send_data(id_flux_oceanq, ocean_qflux, Time)

deallocate(gamma_t            )
deallocate(gamma_q            )
deallocate(en_t               )
deallocate(en_q               )
deallocate(fn_t               )
deallocate(fn_q               )
deallocate(alpha_t            )
deallocate(alpha_q            )
deallocate(alpha_lw           )
deallocate(beta_t             )
deallocate(beta_q             )
deallocate(beta_lw            )
deallocate(eff_heat_capacity  )
deallocate(corrected_flux     )
deallocate(t_surf_dependence  )
!deallocate(delta_t_surf       )

end subroutine mixed_layer

!===================================================================================================

subroutine mixed_layer_end(t_surf,myThid)

real, intent(inout), dimension(:,:) :: t_surf
integer, intent(in)                :: myThid

!integer:: unit

if(.not.module_is_initialized) return

! write a restart file for the surface temperature
!unit = open_file ('RESTART/mixed_layer.res', &
!                  form='native', action='write')
!call write_data (unit, t_surf)
!call close_File (unit)

module_is_initialized = .false.

end subroutine mixed_layer_end

!===================================================================================================

end module mixed_layer_mod