! 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