ecearth3/sources/tm5mp/proj/domino/MOdinHNO3.F90

! First include the set of model-wide compiler flags 
#include "tm5.inc"

! macro
#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN if (status/=0) then; TRACEBACK; return; end if
#define IF_ERROR_RETURN if (status> 0) then; TRACEBACK; return; end if

! First include the set of model-wide compiler flags 
! #include "tm5.inc"

! The ODIN LUT netcdf file contains:
!
! netcdf odin_hno3_o3_ratio {
! dimensions:
! 	nmonth = 12 ;
! 	nlat = 180 ;
! 	npres = 58 ;
! 	nprestm5 = 6 ;
! variables:
! 	float month(nmonth) ;
! 		month:long_name = "index of the month" ;
! 		month:units = "-" ;
! 	float latitude(nlat) ;
! 		latitude:long_name = "latitudes" ;
! 		latitude:units = "degrees north" ;
! 	float pressure(npres) ;
! 		pressure:long_name = "pressure" ;
! 		pressure:units = "hPa" ;
! 	float pressure_tm5(nprestm5) ;
! 		pressure_tm5:long_name = "pressure at the bottom of the TM5 layer" ;
! 		pressure_tm5:units = "hPa" ;
! 	float layerindex_tm5(nprestm5) ;
! 		layerindex_tm5:long_name = "index of the TM5 layer (between 1 and 34)" ;
! 		layerindex_tm5:units = "-" ;
! 	float hno3_o3_ratio(nmonth, nlat, npres) ;
! 		hno3_o3_ratio:long_name = "hno3/o3 ratio" ;
! 		hno3_o3_ratio:units = "ppm/ppm" ;
! 	float hno3_o3_ratio_tm5(nmonth, nlat, nprestm5) ;
! 		hno3_o3_ratio_tm5:long_name = "hno3/o3 ratio integrated over the TM5 layers (34 layer model version)" ;
! 		hno3_o3_ratio_tm5:units = "ppm/ppm" ;
! 	float hno3(nmonth, nlat, npres) ;
! 		hno3:long_name = "mole_fraction_of_hno3_in_air" ;
! 		hno3:units = "mole/mole" ;
! 	float o3(nmonth, nlat, npres) ;
! 		o3:long_name = "mole_fraction_of_ozone_in_air" ;
! 		o3:units = "mole/mole" ;
!
!  pressure_tm5 = 60.18, 39.603, 16.806, 7.132, 2.985, 0.956 ;
!  layerindex_tm5 = 29, 30, 31, 32, 33, 34 ;
!

module MOdinHNO3

  !
  ! Module to load the ODIN climatologies of O3 and HNO3
  !
  !  subroutine ReadOdinHNO3
  !
  !                             Henk Eskes, KNMI, May 2016

  use GO, only : gol, goErr

  implicit none

  private

  public  ::  ReadOdinHNO3, NudgeOdinHNO3, DoneOdinHNO3

  ! Type definition for the ODIN data
  type TOdinHNO3

     integer :: npresTM5    ! nr of TM5 layers on which the climatology is defined
     integer :: nmonth      ! nr of months
     integer :: nlat        ! nr of latitudes

     real, dimension(:), allocatable     :: latitudes          ! the climatology
     real, dimension(:), allocatable     :: pressure_tm5       ! the pressures
     integer, dimension(:), allocatable  :: layerindex_tm5     ! the TM5 layer index
     real, dimension(:,:,:), allocatable :: hno3_o3_ratio_tm5  ! the climatology for TM5 layers

  end type TOdinHNO3

  ! Permanent storage for the ODIN data
  type(TOdinHNO3)  ::  OdinHNO3

  ! Nudging time scales in days
  ! real, dimension(6), parameter  ::  nudgingTime = (/ 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417 /)  ! 1 hour
  ! real, dimension(6), parameter  ::  nudgingTime = (/ 0.25, 0.25, 0.25, 0.25, 0.25, 0.1 /)  ! 6 hours
  real, dimension(6), parameter  ::  nudgingTime = (/ 120.0, 60.0, 30.0, 7.0, 2.0, 1.0 /)
  !                   applied to pressure layers =     60.2, 39.6, 16.8, 7.1, 3.0, 1.0 hPa

  character(len=*), parameter  ::  mname = 'MOdinHNO3Read'

  ! for debug output to screen
  logical, parameter  :: showDebugPrints = .true.

contains



  subroutine NudgeOdinHNO3 ( im, jm, lm, idate, tstep, mixratio_zonal_hno3, mixratio_zonal_o3, rm_hno3 )
    !=======================================================================
    !
    ! Determine scaled HNO3 field based on ODIN climatology
    !
    ! input:
    !   im           = nr of TM5 longitudes
    !   jm           = nr of TM5 latitudes
    !   lm           = nr of TM5 layers (climatol only defined for lm=34)
    !   idate        = current date-time
    !   tstep        = time step in seconds (integer)
    !   mixratio_zonal_hno3 = zonal-mean hno3 mixing ratio in TM5
    !   mixratio_zonal_o3   = zonal-mean ozone mixing ratio in TM5
    ! input/output:
    !   rm_hno3      = HNO3 mass in TM5 
    !
    ! The routine will allocate the arrays in "OdinHNO3",
    !     using the dimensions in the file  
    !                                            Henk Eskes, KNMI, May 2016
    !=======================================================================

    implicit none

    ! in/out
    integer, intent(in)                       ::  im, jm, lm
    integer, intent(in)                       ::  tstep
    integer, dimension(6), intent(in)         ::  idate
    real, dimension(jm,lm), intent(in)        ::  mixratio_zonal_hno3
    real, dimension(jm,lm), intent(in)        ::  mixratio_zonal_o3
    real, dimension(im,jm,lm), intent(inout)  ::  rm_hno3

    ! local
    integer  ::  i, j, l, lOdin
    real     ::  c, tfac, expfac, ratio_model, ratio_odin_intp, ratio
    integer  ::  month, month_other
    real     ::  fraction_othermonth, fraction_thismonth

    ! start code

    ! determine time interpolation factors between 2 months
    month = idate(2) 
    month_other = month
    fraction_thismonth = 1.0
    fraction_othermonth = 0.0 
    if ( idate(3) .lt. 10 ) then
       month_other = month - 1
       if ( month_other < 1 ) month_other = 12
       fraction_othermonth = ( 10.0 - real(idate(3)) ) / 19.0
       fraction_thismonth = 1.0 - fraction_othermonth
    end if
    if ( idate(3) .gt. 20 ) then
       month_other = month + 1
       if ( month_other > 12 ) month_other = 1
       fraction_othermonth = (real(idate(3))-20.0) / 22.0
       fraction_thismonth = 1.0 - fraction_othermonth
    end if

    ! nudge HNO3 field to ratio of ODIN HNO3/O3 climatology using the zonal mean
    do lOdin = 1, OdinHNO3%npresTM5
       tfac = real(tstep) / (nudgingTime(lOdin)*3600.0*24.0)   ! nudgingTime in days, tstep in seconds
       expfac = exp(-1.0*tfac)
       do j = 1, jm
          l = OdinHNO3%layerindex_tm5(lOdin)
          ratio_odin_intp = fraction_thismonth  * OdinHNO3%hno3_o3_ratio_tm5(l,j,month) + &
                            fraction_othermonth * OdinHNO3%hno3_o3_ratio_tm5(l,j,month_other) 
          ratio_model = mixratio_zonal_hno3(j,l) / mixratio_zonal_o3(j,l) 
          ratio = ratio_odin_intp/ratio_model
          c = ratio - expfac*(ratio_odin_intp-ratio_model)/ratio_model
          if ( showDebugPrints ) then
             if ( j == 120 ) then
                if ( lOdin == 1 ) print '(a,i5,5i3)', 'idate = ',idate
                print *, ''
                print '(a,5i4,f8.3)', 'j, l, lOdin, month, month_other, pres        = ',j, l, lOdin, month, month_other, OdinHNO3%pressure_tm5(lOdin)
                print '(a,4e12.3)',   'Odin: intp, month, other, fraction_thismonth = ',ratio_odin_intp,Odinhno3%hno3_o3_ratio_tm5(l,j,month),OdinHNO3%hno3_o3_ratio_tm5(l,j,month_other), fraction_thismonth
                print '(a,3e12.3)',   'hno3/o3, hno3, o3                            = ', mixratio_zonal_hno3(j,l)/mixratio_zonal_o3(j,l) , mixratio_zonal_hno3(j,l), mixratio_zonal_o3(j,l)
                print '(a,f10.5,f7.0,f10.6,f7.1,f10.6)', 'ratio, tstep, tfac, nudgingTime, multfactor  = ',ratio, real(tstep), tfac, nudgingTime(lOdin), c
             end if
          end if
          ! multiply all zonal values with the same correction factor
          do i = 1, im
             rm_hno3(i,j,l) = rm_hno3(i,j,l) * c
          end do
       end do
    end do
    if ( showDebugPrints ) then
       print *,   ' '
       print '(a,6f8.3)', 'nudging times = ',nudgingTime 
       print '(a,4e12.3)', 'min/max of hno3, o3                          = ', minval(mixratio_zonal_hno3(:,:)), maxval(mixratio_zonal_hno3(:,:)), minval(mixratio_zonal_o3(:,:)), maxval(mixratio_zonal_o3(:,:))
    end if

  end subroutine NudgeOdinHNO3



  subroutine ReadOdinHNO3 ( OdinHNO3File, jm, lm, status )
    !=======================================================================
    !
    ! ReadOdinHNO3:  read ODIN climatology of [HNO3]/[O3]
    !                as a function of the month and the TM5 layer
    !                (lookup table created by Henk Eskes, KNMI, April 2016)
    !
    ! input:
    !   OdinHNO3File = path+name of the file containing the ODIN Climatology
    !   jm           = nr of TM5 latitudes
    !   lm           = nr of TM5 layers (climatol only defined for lm=34)
    ! output:
    !   status = status ( 0 means OK )
    !
    ! The routine will allocate the arrays in "OdinHNO3",
    !     using the dimensions in the file  
    !                                            Henk Eskes, KNMI, May 2016
    !=======================================================================

    use MDF,      only : MDF_Open,      MDF_NETCDF,  MDF_READ 
    use MDF,      only : MDF_Inq_VarID, MDF_Get_Var, MDF_Close 
    use MDF,      only : MDF_Inquire_Dimension, MDF_Inq_DimID

    implicit none

    ! in/out
    character(*), intent(in)        :: OdinHNO3File
    integer, intent(in)             :: jm, lm
    integer, intent(out)            :: status

    ! local
    character(len=*), parameter ::  rname = mname//'/ReadOdinHNO3'

    integer      :: fid
    integer      :: dimid, id_lat, id_pres, id_ind, id_hno3

    integer, dimension(:), allocatable :: dphi_int, psurf_int

    real :: amfGeo

    ! start code

    status = 0

    call MDF_Open( trim(OdinHNO3File), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN

    ! get the dimensions
    call MDF_Inq_DimID( fid, 'nmonth', dimid, status )
    IF_NOTOK_RETURN
    call MDF_Inquire_Dimension( fid, dimid, status, length=OdinHNO3%Nmonth )
    IF_NOTOK_RETURN
    call MDF_Inq_DimID( fid, 'nlat', dimid, status )
    IF_NOTOK_RETURN
    call MDF_Inquire_Dimension( fid, dimid, status, length=OdinHNO3%nlat )
    IF_NOTOK_RETURN
    call MDF_Inq_DimID( fid, 'nprestm5', dimid, status )
    IF_NOTOK_RETURN
    call MDF_Inquire_Dimension( fid, dimid, status, length=OdinHNO3%nprestm5 )
    IF_NOTOK_RETURN

    print '(a,3i4)', 'ReadOdinHNO3, array dimensions (pressure,lat,month) = ', OdinHNO3%nprestm5, OdinHNO3%nlat, OdinHNO3%nmonth
       
    ! Allocate storage space
    allocate ( OdinHNO3%latitudes(OdinHNO3%nlat) )
    allocate ( OdinHNO3%pressure_tm5(OdinHNO3%nprestm5) )
    allocate ( OdinHNO3%layerindex_tm5(OdinHNO3%nprestm5) )
    allocate ( OdinHNO3%hno3_o3_ratio_tm5(OdinHNO3%nprestm5, OdinHNO3%nlat, OdinHNO3%nmonth) )

    ! Define variables
    CALL MDF_Inq_VarID( fid, 'latitude', id_lat, status )
    IF_ERROR_RETURN
    CALL MDF_Inq_VarID( fid, 'pressure_tm5', id_pres, status )
    IF_ERROR_RETURN
    CALL MDF_Inq_VarID( fid, 'layerindex_tm5', id_ind, status )
    IF_ERROR_RETURN
    CALL MDF_Inq_VarID( fid, 'hno3_o3_ratio_tm5', id_hno3, status )
    IF_ERROR_RETURN

    ! Read variables
    CALL MDF_Get_Var( fid, id_lat, OdinHNO3%latitudes, status )      ! latitudes
    IF_NOTOK_RETURN

    CALL MDF_Get_Var( fid, id_pres, OdinHNO3%pressure_tm5, status )    ! pressures
    IF_NOTOK_RETURN

    CALL MDF_Get_Var( fid, id_ind, OdinHNO3%layerindex_tm5, status )  ! TM5 layer indices
    IF_NOTOK_RETURN

    CALL MDF_Get_Var( fid, id_hno3, OdinHNO3%hno3_o3_ratio_tm5, status )     ! HNO3 / O3 ratios at the TM5 layers
    IF_NOTOK_RETURN

    ! Close file
    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN

    ! check if this version of TM5 is consistent with the climatology
    if ( (jm /= OdinHNO3%nlat) .or. (lm /= OdinHNO3%layerindex_tm5(OdinHNO3%nprestm5)) ) then
       status = -1
       print *, 'ReadOdinHNO3: ERROR dimensions incorrect'
       print *, '   TM5 jm, lm = ',jm,lm
       print *, '   ODIN climatol jm, lm = ',OdinHNO3%nlat,OdinHNO3%layerindex_tm5(OdinHNO3%nprestm5)
       IF_NOTOK_RETURN
    end if

    if ( showDebugPrints ) then
       print '(a,2f8.2)', 'ReadOdinHNO3, lat range  = ',OdinHNO3%latitudes(1),OdinHNO3%latitudes(OdinHNO3%nlat)
       print '(a,6f10.4)', 'ReadOdinHNO3, press = ',OdinHNO3%pressure_tm5
       print '(a,6i4)', 'ReadOdinHNO3, index = ',OdinHNO3%layerindex_tm5
       print '(a,6e12.3)', 'ReadOdinHNO3, hno3/o3, month 12, lat 120 = ',OdinHNO3%hno3_o3_ratio_tm5(:,120,12)
    end if

  end subroutine ReadOdinHNO3


  subroutine DoneOdinHNO3 ( )
    !
    ! Deallocate arrays
    !

    implicit none

    ! De-allocate storage space of the ODIN climatology
    if ( allocated ( OdinHNO3%latitudes ) ) deallocate ( OdinHNO3%latitudes )
    if ( allocated ( OdinHNO3%pressure_tm5 ) ) deallocate ( OdinHNO3%pressure_tm5 )
    if ( allocated ( OdinHNO3%layerindex_tm5 ) ) deallocate ( OdinHNO3%layerindex_tm5 )
    if ( allocated ( OdinHNO3%hno3_o3_ratio_tm5 ) ) deallocate ( OdinHNO3%hno3_o3_ratio_tm5 )

  end subroutine DoneOdinHNO3


end module MOdinHNO3