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

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, c20, &
+      use ice_constants, only: c0, c1, c2, c3, c4, c5, c8, c10, c12, c20, &
                                c180, c365, c1000, c3600
       use ice_constants, only: p001, p01, p1, 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'
@@ -244,6 +245,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 +542,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 +1381,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_format) == 'nc') 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 +2011,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 +2250,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 6 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.jra55

@@ -0,0 +1,20 @@
+dt = 3600.0
+runtype = 'initial'
+ice_ic  = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc'
+grid_format = 'bin'
+grid_type   = 'displaced_pole'
+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'
+maskhalo_dyn   = .true.
+maskhalo_remap = .true.
+maskhalo_bound = .true.
+fyear_init      = 2005
+ycycle          = 1
+atm_data_dir    = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55'
+atm_data_format = 'nc'
+atm_data_type   = 'JRA55'
+precip_units    = 'mm_per_sec'
+ocn_data_dir    = 'default'
+bgc_data_dir    = 'default'
+

configuration/scripts/tests/base_suite.ts

@@ -36,4 +36,4 @@ smoke          gx3     8x2        bgcz,debug
 smoke          gx3     8x1        bgcskl,debug
 #smoke          gx3     4x1        bgcz,thread        smoke_gx3_8x2_bgcz
 restart        gx1     4x2        bgcsklclim,medium
-restart        gx1     8x1        bgczclim,medium
+restart        gx1     8x1        jra55