Added new capability in ice_forcing.F90 to test with JRA55 atmospheri… (

#350)

* Added new capability in ice_forcing.F90 to test with JRA55 atmospheric forcing on gx1 grid.

* Added JRA55 as an option under atm_data_type

* Modified ice_forcing.F90 by replacing
elseif (trim(atm_data_format) == 'nc') then ! netcdf

with

elseif (trim(atm_data_type) == 'JRA55') then ! netcdf

Corrected comment regarding JRA55 forcing to state "E.g. record 1 gives conditions at 3 am GMT on 1 January"

Modified ./scripts/options/set_nml.jra55 by
setting use_restart_time = .false.
setting year_init = 2005
setting fyear_init - 2005
and corrected units for precip_units    = 'mks'

edited ./scripts/tests/base_suite.ts by reinstating bgc restart test and added a jra55 "short" smoke test.

* fix an incorrect test in the base suite

* Add a check in init_forcing_atmo to check if 'use_leap_years' is set to true.  If so, abort if the atm_data_type is not JRA55, default, or box2001

* Fix a mistake in diagnostic warning output for the use_leap_years check

* set use_leap_years to false in the boxrestore namelist option file

* add use_restart_time specification to test_nml.restart* files

* Remove options in JRA55 namelist option file that are the same as the default in ice_in.  Also remove the specification of use_restart_time to .false.

* fix a bug in the newly-implemented use_leap_years check during initialization

* add use_leap_years and use_restart_time specification to jra55 namelist option file

* add new jra55_2008 and run90day namelist option files

* modify the JRA55 smoke test to be a 90 day test and start in 2008.  Also add a restart test for JRA55 in base_suite

* update documentation to include citation for JRA55

cicecore/cicedynB/general/ice_forcing.F90

@@ -30,7 +30,7 @@ module ice_forcing
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite, &
                             timer_bound
       use ice_arrays_column, only: oceanmixed_ice, restore_bgc
-      use ice_constants, only: c0, c1, c2, c3, c4, c5, c10, c12, c15, c20, &
+      use ice_constants, only: c0, c1, c2, c3, c4, c5, c8, c10, c12, c15, c20, &
                                c180, c360, c365, c1000, c3600
       use ice_constants, only: p001, p01, p1, p2, p25, p5, p6
       use ice_constants, only: cm_to_m
@@ -115,7 +115,8 @@ module ice_forcing
          atm_data_format, & ! 'bin'=binary or 'nc'=netcdf
          ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf
          atm_data_type, & ! 'default', 'monthly', 'ncar', 
-                          ! 'LYq' or 'hadgem' or 'oned'
+                          ! 'LYq' or 'hadgem' or 'oned' or
+                          ! 'JRA55'
          bgc_data_type, & ! 'default', 'clim'
          ocn_data_type, & ! 'default', 'clim', 'ncar', 'oned',
                           ! 'hadgem_sst' or 'hadgem_sst_uvocn'
@@ -211,6 +212,8 @@ subroutine init_forcing_atmo
 ! Determine the current and final year of the forcing cycle based on
 ! namelist input; initialize the atmospheric forcing data filenames.
 
+      use ice_calendar, only: use_leap_years
+
       character(len=*), parameter :: subname = '(init_forcing_atmo)'
 
       ! Allocate forcing arrays 
@@ -235,6 +238,14 @@ subroutine init_forcing_atmo
             file=__FILE__, line=__LINE__)
       endif
 
+      if (use_leap_years .and. (trim(atm_data_type) /= 'JRA55' .and. &
+                                trim(atm_data_type) /= 'default' .and. &
+                                trim(atm_data_type) /= 'box2001')) then
+         write(nu_diag,*) 'use_leap_years option is currently only supported for'
+         write(nu_diag,*) 'JRA55, default , and box2001 atmospheric data'
+         call abort_ice(error_message=subname, file=__FILE__, line=__LINE__)
+      endif
+
     !-------------------------------------------------------------------
     ! Get filenames for input forcing data     
     !-------------------------------------------------------------------
@@ -244,6 +255,8 @@ subroutine init_forcing_atmo
          call NCAR_files(fyear)
       elseif (trim(atm_data_type) == 'LYq') then
          call LY_files(fyear)
+      elseif (trim(atm_data_type) == 'JRA55') then
+         call JRA55_files(fyear)
       elseif (trim(atm_data_type) == 'hadgem') then
          call hadgem_files(fyear)
       elseif (trim(atm_data_type) == 'monthly') then
@@ -539,6 +552,8 @@ subroutine get_forcing_atmo
          call ncar_data
       elseif (trim(atm_data_type) == 'LYq') then
          call LY_data
+      elseif (trim(atm_data_type) == 'JRA55') then
+         call JRA55_data(fyear)
       elseif (trim(atm_data_type) == 'hadgem') then
          call hadgem_data
       elseif (trim(atm_data_type) == 'monthly') then
@@ -1376,6 +1391,10 @@ subroutine file_year (data_file, yr)
          i = index(data_file,'.nc') - 5
          tmpname = data_file
          write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc'
+      elseif (trim(atm_data_type) == 'JRA55') then ! netcdf
+         i = index(data_file,'.nc') - 5
+         tmpname = data_file
+         write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc'
       else                                     ! LANL/NCAR naming convention
          i = index(data_file,'.dat') - 5
          tmpname = data_file
@@ -2002,6 +2021,22 @@ subroutine LY_files (yr)
       endif                     ! master_task
 
       end subroutine LY_files
+      subroutine JRA55_files(yr)
+!
+      integer (kind=int_kind), intent(in) :: &
+           yr                   ! current forcing year
+
+      character(len=*), parameter :: subname = '(JRA55_files)'
+
+      uwind_file = &
+           trim(atm_data_dir)//'/8XDAILY/JRA55_03hr_forcing_2005.nc'
+      call file_year(uwind_file,yr)
+  if (my_task == master_task) then
+         write (nu_diag,*) ' '
+         write (nu_diag,*) 'Atmospheric data files:'
+         write (nu_diag,*) trim(uwind_file)
+    endif
+      end subroutine JRA55_files
 
 !=======================================================================
 !
@@ -2225,6 +2260,239 @@ subroutine LY_data
 
       end subroutine LY_data
 
+      subroutine JRA55_data (yr)
+
+      use ice_blocks, only: block, get_block
+      use ice_global_reductions, only: global_minval, global_maxval
+      use ice_domain, only: nblocks, distrb_info, blocks_ice
+      use ice_flux, only: fsnow, Tair, uatm, vatm, Qa, fsw, flw
+      use ice_grid, only: hm, tlon, tlat, tmask, umask
+      use ice_state, only: aice
+      use ice_calendar, only: days_per_year, use_leap_years
+
+
+      integer (kind=int_kind) :: & 
+          ncid        , & ! netcdf file id
+	  i, j        , &
+          ixm,ixx,ixp , & ! record numbers for neighboring months
+          recnum      , & ! record number
+          maxrec      , & ! maximum record number
+          recslot     , & ! spline slot for current record
+          midmonth    , & ! middle day of month
+          dataloc     , & ! = 1 for data located in middle of time interval
+                          ! = 2 for date located at end of time interval
+          iblk        , & ! block index
+          ilo,ihi,jlo,jhi, & ! beginning and end of physical domain
+          yr              ! current forcing year
+
+      real (kind=dbl_kind) :: &
+          sec3hr          , & ! number of seconds in 3 hours
+          secday          , & ! number of seconds in day
+          Tffresh         , &
+          vmin, vmax
+
+      logical (kind=log_kind) :: readm, read6,debug_n_d
+
+      type (block) :: &
+         this_block           ! block information for current block
+
+      character(len=64) :: fieldname !netcdf field name
+      character(len=*), parameter :: subname = '(JRA55_data)'
+
+      debug_n_d = .false.  !usually false
+
+      call icepack_query_parameters(Tffresh_out=Tffresh)
+      call icepack_query_parameters(secday_out=secday)
+      call icepack_warnings_flush(nu_diag)
+      if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+         file=__FILE__, line=__LINE__)
+
+    !-------------------------------------------------------------------
+    ! 3-hourly data
+    ! 
+    ! Assume that the 3-hourly value is located at the end of the
+    !  3-hour period.  This is the convention for NCEP reanalysis data.
+    !  E.g. record 1 gives conditions at 3 am GMT on 1 January.
+    !-------------------------------------------------------------------
+
+      dataloc = 2               ! data located at end of interval
+      sec3hr = secday/c8        ! seconds in 3 hours
+      !maxrec = 2920            ! 365*8; for leap years = 366*8
+
+      if(use_leap_years) days_per_year = 366 !overrides setting of 365 in ice_calendar
+      maxrec = days_per_year*8
+
+      if(days_per_year == 365 .and. (mod(yr,  4) == 0)) then
+      call abort_ice('days_per_year should be set to 366 for leap years')
+      end if
+
+      ! current record number
+      recnum = 8*int(yday) - 7 + int(real(sec,kind=dbl_kind)/sec3hr)
+
+      ! Compute record numbers for surrounding data (2 on each side)
+
+      ixm = mod(recnum+maxrec-2,maxrec) + 1
+      ixx = mod(recnum-1,       maxrec) + 1
+
+      ! Compute interpolation coefficients
+      ! If data is located at the end of the time interval, then the
+      !  data value for the current record goes in slot 2
+
+      recslot = 2
+      ixp = -99
+      call interp_coeff (recnum, recslot, sec3hr, dataloc)
+
+      ! Read
+      read6 = .false.
+      if (istep==1 .or. oldrecnum .ne. recnum) read6 = .true.
+         !-------------------------------------------------------------------
+         ! File is NETCDF with winds in NORTH and EAST direction
+         ! file variable names are:
+         ! glbrad   (shortwave W/m^2)
+         ! dlwsfc   (longwave W/m^2)
+         ! wndewd   (eastward wind m/s)
+         ! wndnwd   (northward wind m/s)
+         ! airtmp   (air temperature K)
+         ! spchmd   (specific humidity kg/kg)
+         ! ttlpcp   (precipitation kg/m s-1)
+         !-------------------------------------------------------------------
+         call ice_open_nc(uwind_file,ncid)
+
+	 fieldname = 'airtmp'
+         call ice_read_nc(ncid,recnum,fieldname,Tair_data(:,:,1,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+         call ice_read_nc(ncid,recnum,fieldname,Tair_data(:,:,2,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+
+	 fieldname = 'wndewd'
+         call ice_read_nc(ncid,recnum,fieldname,uatm_data(:,:,1,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+         call ice_read_nc(ncid,recnum,fieldname,uatm_data(:,:,2,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+
+	 fieldname = 'wndnwd'
+         call ice_read_nc(ncid,recnum,fieldname,vatm_data(:,:,1,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+         call ice_read_nc(ncid,recnum,fieldname,vatm_data(:,:,2,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+
+	 fieldname = 'spchmd'
+         call ice_read_nc(ncid,recnum,fieldname,Qa_data(:,:,1,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+         call ice_read_nc(ncid,recnum,fieldname,Qa_data(:,:,2,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+
+	 fieldname = 'glbrad'
+         call ice_read_nc(ncid,recnum,fieldname,fsw_data(:,:,1,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+         call ice_read_nc(ncid,recnum,fieldname,fsw_data(:,:,2,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+
+	 fieldname = 'dlwsfc'
+         call ice_read_nc(ncid,recnum,fieldname,flw_data(:,:,1,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+         call ice_read_nc(ncid,recnum,fieldname,flw_data(:,:,2,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+
+	 fieldname = 'ttlpcp'
+         call ice_read_nc(ncid,recnum,fieldname,fsnow_data(:,:,1,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+         call ice_read_nc(ncid,recnum,fieldname,fsnow_data(:,:,2,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+
+	 call ice_close_nc(ncid)
+
+
+      ! Interpolate
+      call interpolate_data (Tair_data, Tair)
+      call interpolate_data (uatm_data, uatm)
+      call interpolate_data (vatm_data, vatm)
+      call interpolate_data (Qa_data, Qa)
+      call interpolate_data (fsw_data, fsw)
+      call interpolate_data (flw_data, flw)
+      call interpolate_data (fsnow_data, fsnow)
+
+      !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
+      do iblk = 1, nblocks
+        ! limit summer Tair values where ice is present
+        do j = 1, ny_block
+          do i = 1, nx_block
+            if (aice(i,j,iblk) > p1) Tair(i,j,iblk) = min(Tair(i,j,iblk), Tffresh+p1)
+          enddo
+        enddo
+
+        do j = 1, ny_block
+          do i = 1, nx_block
+            Qa  (i,j,iblk) = Qa  (i,j,iblk) * hm(i,j,iblk)
+            Tair(i,j,iblk) = Tair(i,j,iblk) * hm(i,j,iblk)
+            uatm(i,j,iblk) = uatm(i,j,iblk) * hm(i,j,iblk)
+            vatm(i,j,iblk) = vatm(i,j,iblk) * hm(i,j,iblk)
+            fsw (i,j,iblk) = fsw (i,j,iblk) * hm(i,j,iblk)
+            flw (i,j,iblk) = flw (i,j,iblk) * hm(i,j,iblk)
+            fsnow(i,j,iblk) = fsnow (i,j,iblk) * hm(i,j,iblk)
+          enddo
+        enddo
+
+      enddo  ! iblk
+      !$OMP END PARALLEL DO
+
+      ! Save record number
+      oldrecnum = recnum
+
+         if (dbug) then
+           if (my_task == master_task) write (nu_diag,*) 'JRA55_bulk_data'
+           vmin = global_minval(fsw,distrb_info,tmask)
+
+           vmax = global_maxval(fsw,distrb_info,tmask)
+           if (my_task.eq.master_task)  &
+               write (nu_diag,*) 'fsw',vmin,vmax 
+           vmin = global_minval(flw,distrb_info,tmask)
+           vmax = global_maxval(flw,distrb_info,tmask)
+           if (my_task.eq.master_task) & 
+               write (nu_diag,*) 'flw',vmin,vmax
+           vmin =global_minval(fsnow,distrb_info,tmask)
+           vmax =global_maxval(fsnow,distrb_info,tmask)
+           if (my_task.eq.master_task) & 
+               write (nu_diag,*) 'fsnow',vmin,vmax
+           vmin = global_minval(Tair,distrb_info,tmask)
+           vmax = global_maxval(Tair,distrb_info,tmask)
+           if (my_task.eq.master_task) & 
+               write (nu_diag,*) 'Tair',vmin,vmax
+           vmin = global_minval(uatm,distrb_info,umask)
+           vmax = global_maxval(uatm,distrb_info,umask)
+           if (my_task.eq.master_task) & 
+               write (nu_diag,*) 'uatm',vmin,vmax
+           vmin = global_minval(vatm,distrb_info,umask)
+           vmax = global_maxval(vatm,distrb_info,umask)
+           if (my_task.eq.master_task) & 
+               write (nu_diag,*) 'vatm',vmin,vmax
+           vmin = global_minval(Qa,distrb_info,tmask)
+           vmax = global_maxval(Qa,distrb_info,tmask)
+           if (my_task.eq.master_task)  &
+               write (nu_diag,*) 'Qa',vmin,vmax
+           if (my_task.eq.master_task)  &
+               write (nu_diag,*) 'maxrec',maxrec
+               write (nu_diag,*) 'days_per_year', days_per_year
+
+        endif                   ! dbug
+
+      end subroutine JRA55_data
+
+
 !=======================================================================
 !
 ! AOMIP shortwave forcing

configuration/scripts/options/set_nml.boxrestore

@@ -2,7 +2,7 @@ nilyr =	1
 ice_ic = 'default'
 restart = .false.
 restart_ext = .true.
-use_leap_years = .true.
+use_leap_years = .false.
 ndtd = 2
 kcatbound = 1
 ew_boundary_type = 'cyclic'

configuration/scripts/options/set_nml.jra55

@@ -0,0 +1,17 @@
+year_init      = 2005
+ice_ic  = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc'
+grid_file   = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin'
+kmt_file    = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin'
+bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc'
+use_leap_years = .true.
+use_restart_time = .false.
+maskhalo_dyn   = .true.
+maskhalo_remap = .true.
+maskhalo_bound = .true.
+fyear_init      = 2005
+atm_data_dir    = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55'
+atm_data_format = 'nc'
+atm_data_type   = 'JRA55'
+precip_units    = 'mks'
+ocn_data_dir    = 'default'
+bgc_data_dir    = 'default'

configuration/scripts/options/set_nml.jra55_2008

@@ -0,0 +1,17 @@
+year_init      = 2008
+ice_ic  = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc'
+grid_file   = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin'
+kmt_file    = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin'
+bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc'
+use_leap_years = .true.
+use_restart_time = .false.
+maskhalo_dyn   = .true.
+maskhalo_remap = .true.
+maskhalo_bound = .true.
+fyear_init      = 2008
+atm_data_dir    = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55'
+atm_data_format = 'nc'
+atm_data_type   = 'JRA55'
+precip_units    = 'mks'
+ocn_data_dir    = 'default'
+bgc_data_dir    = 'default'

configuration/scripts/options/set_nml.run90day

@@ -0,0 +1,5 @@
+npt            = 2160
+dumpfreq       = 'd'
+dumpfreq_n     = 30
+histfreq       = 'd','x','x','x','x'
+f_aice	       = 'd'

configuration/scripts/options/test_nml.restart1

@@ -2,3 +2,4 @@ npt            = 240
 dumpfreq       = 'd'
 dumpfreq_n     = 5
 runtype        = 'initial'
+use_restart_time = .false.

configuration/scripts/options/test_nml.restart2

@@ -3,4 +3,4 @@ dumpfreq       = 'd'
 dumpfreq_n     = 5
 runtype        = 'continue'
 restart        = .true.
-
+use_restart_time = .true.