Merge pull request #2 from CICE-Consortium/master

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

@@ -85,8 +85,14 @@ module ice_dyn_shared
       logical (kind=log_kind), public :: &
          basalstress   ! if true, basal stress for landfast on
 
+      ! basal stress parameters
       real (kind=dbl_kind), public :: &
-         k1            ! 1st free parameter for landfast parameterization
+         k1, &        ! 1st free parameter for landfast parameterization
+         k2, &        ! second free parameter (N/m^3) for landfast parametrization 
+         alphab, &    ! alphab=Cb factor in Lemieux et al 2015
+         threshold_hw ! max water depth for grounding 
+                      ! see keel data from Amundrud et al. 2004 (JGR)
+
 
 !=======================================================================
 
@@ -887,11 +893,7 @@ subroutine basal_stress_coeff (nx_block, ny_block,         &
          au,  & ! concentration of ice at u location
          hu,  & ! volume per unit area of ice at u location (mean thickness)
          hwu, & ! water depth at u location
-         hcu, & ! critical thickness at u location
-         k2 = 15.0_dbl_kind , &     ! second free parameter (N/m^3) for landfast parametrization 
-         alphab = 20.0_dbl_kind, &  ! alphab=Cb factor in Lemieux et al 2015
-         threshold_hw = 30.0_dbl_kind ! max water depth for grounding 
-                                      ! see keel data from Amundrud et al. 2004 (JGR)
+         hcu    ! critical thickness at u location
 
       integer (kind=int_kind) :: &
          i, j, ij

cicecore/cicedynB/general/ice_forcing.F90

@@ -117,7 +117,7 @@ module ice_forcing
          ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf
          atm_data_type, & ! 'default', 'monthly', 'ncar', 
                           ! 'LYq' or 'hadgem' or 'oned' or
-                          ! 'JRA55'
+                          ! 'JRA55_gx1' or 'JRA55_gx3'
          bgc_data_type, & ! 'default', 'clim'
          ocn_data_type, & ! 'default', 'clim', 'ncar', 'oned',
                           ! 'hadgem_sst' or 'hadgem_sst_uvocn'
@@ -241,8 +241,8 @@ 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. &
+      if (use_leap_years .and. (trim(atm_data_type) /= 'JRA55_gx1' .and. &
+                                trim(atm_data_type) /= 'JRA55_gx3' .and. &
                                 trim(atm_data_type) /= 'hycom' .and. &
                                 trim(atm_data_type) /= 'box2001')) then
          write(nu_diag,*) 'use_leap_years option is currently only supported for'
@@ -259,8 +259,10 @@ 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) == 'JRA55_gx1') then
+         call JRA55_gx1_files(fyear)
+      elseif (trim(atm_data_type) == 'JRA55_gx3') then
+         call JRA55_gx3_files(fyear)
       elseif (trim(atm_data_type) == 'hadgem') then
          call hadgem_files(fyear)
       elseif (trim(atm_data_type) == 'monthly') then
@@ -556,7 +558,9 @@ subroutine get_forcing_atmo
          call ncar_data
       elseif (trim(atm_data_type) == 'LYq') then
          call LY_data
-      elseif (trim(atm_data_type) == 'JRA55') then
+      elseif (trim(atm_data_type) == 'JRA55_gx1') then
+         call JRA55_data(fyear)
+      elseif (trim(atm_data_type) == 'JRA55_gx3') then
          call JRA55_data(fyear)
       elseif (trim(atm_data_type) == 'hadgem') then
          call hadgem_data
@@ -1395,7 +1399,11 @@ 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
+      elseif (trim(atm_data_type) == 'JRA55_gx1') then ! netcdf
+         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_gx3') then ! netcdf
          i = index(data_file,'.nc') - 5
          tmpname = data_file
          write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc'
@@ -2025,12 +2033,12 @@ subroutine LY_files (yr)
       endif                     ! master_task
 
       end subroutine LY_files
-      subroutine JRA55_files(yr)
+      subroutine JRA55_gx1_files(yr)
 !
       integer (kind=int_kind), intent(in) :: &
            yr                   ! current forcing year
 
-      character(len=*), parameter :: subname = '(JRA55_files)'
+      character(len=*), parameter :: subname = '(JRA55_gx1_files)'
 
       uwind_file = &
            trim(atm_data_dir)//'/8XDAILY/JRA55_03hr_forcing_2005.nc'
@@ -2040,8 +2048,23 @@ subroutine JRA55_files(yr)
          write (nu_diag,*) 'Atmospheric data files:'
          write (nu_diag,*) trim(uwind_file)
     endif
-      end subroutine JRA55_files
+      end subroutine JRA55_gx1_files
+      subroutine JRA55_gx3_files(yr)
+!
+      integer (kind=int_kind), intent(in) :: &
+           yr                   ! current forcing year
+
+      character(len=*), parameter :: subname = '(JRA55_gx3_files)'
 
+      uwind_file = &
+           trim(atm_data_dir)//'/8XDAILY/JRA55_gx3_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_gx3_files
 !=======================================================================
 !
 ! read Large and Yeager atmospheric data

cicecore/cicedynB/general/ice_init.F90

@@ -95,7 +95,8 @@ subroutine input_data
                           dxrect, dyrect
       use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
                                 kevp_kernel, &
-                                basalstress, k1, Ktens, e_ratio, coriolis, &
+                                basalstress, k1, k2, alphab, threshold_hw, &
+                                Ktens, e_ratio, coriolis, &
                                 kridge, ktransport, brlx, arlx
       use ice_transport_driver, only: advection
       use ice_restoring, only: restore_ice
@@ -184,7 +185,7 @@ subroutine input_data
         advection,      coriolis,       kridge,         ktransport,     &
         kstrength,      krdg_partic,    krdg_redist,    mu_rdg,         &
         e_ratio,        Ktens,          Cf,             basalstress,    &
-        k1
+        k1,             k2,             alphab,         threshold_hw      
 
       namelist /shortwave_nml/ &
         shortwave,      albedo_type,                                    &
@@ -296,6 +297,9 @@ subroutine input_data
       close_boundaries = .false.   ! true = set land on edges of grid
       basalstress= .false.   ! if true, basal stress for landfast is on
       k1 = 8.0_dbl_kind      ! 1st free parameter for landfast parameterization
+      k2 = 15.0_dbl_kind     ! dah: second free parameter (N/m^3) for landfast parametrization
+      alphab = 20.0_dbl_kind       ! alphab=Cb factor in Lemieux et al 2015
+      threshold_hw = 30.0_dbl_kind ! max water depth for grounding
       Ktens = 0.0_dbl_kind   ! T=Ktens*P (tensile strength: see Konig and Holland, 2010)
       e_ratio = 2.0_dbl_kind ! EVP ellipse aspect ratio
       advection  = 'remap'   ! incremental remapping transport scheme
@@ -569,6 +573,9 @@ subroutine input_data
       call broadcast_scalar(Cf,                 master_task)
       call broadcast_scalar(basalstress,        master_task)
       call broadcast_scalar(k1,                 master_task)
+      call broadcast_scalar(k2,                 master_task)
+      call broadcast_scalar(alphab,             master_task)
+      call broadcast_scalar(threshold_hw,       master_task)
       call broadcast_scalar(Ktens,              master_task)
       call broadcast_scalar(e_ratio,            master_task)
       call broadcast_scalar(advection,          master_task)
@@ -1036,8 +1043,12 @@ subroutine input_data
          write(nu_diag,1000) ' mu_rdg                    = ', mu_rdg
          if (kstrength == 1) &
          write(nu_diag,1000) ' Cf                        = ', Cf
+
          write(nu_diag,1010) ' basalstress               = ', basalstress
          write(nu_diag,1005) ' k1                        = ', k1
+         write(nu_diag,1005) ' k2                        = ', k2
+         write(nu_diag,1005) ' alphab                    = ', alphab
+         write(nu_diag,1005) ' threshold_hw              = ', threshold_hw
          write(nu_diag,1005) ' Ktens                     = ', Ktens
          write(nu_diag,1005) ' e_ratio                   = ', e_ratio    
          write(nu_diag,1030) ' advection                 = ', &

cicecore/cicedynB/infrastructure/ice_read_write.F90

@@ -28,6 +28,11 @@ module ice_read_write
       implicit none
 
       private
+
+      integer (kind=int_kind), parameter, private :: &
+           bits_per_byte = 8 ! number of bits per byte. 
+                             ! used to determine RecSize in ice_open
+
       public :: ice_open,           &
                 ice_open_ext,       &
                 ice_open_nc,        &
@@ -86,7 +91,7 @@ subroutine ice_open(nu, filename, nbits, algn)
            nbits         ! no. of bits per variable (0 for sequential access)
 
       integer (kind=int_kind), intent(in), optional :: algn
-      integer (kind=int_kind) :: RecSize, Remnant
+      integer (kind=int_kind) :: RecSize, Remnant, nbytes
 
       character (*) :: filename
 
@@ -99,7 +104,13 @@ subroutine ice_open(nu, filename, nbits, algn)
             open(nu,file=filename,form='unformatted')
 
          else                   ! direct access
-            RecSize = nx_global*ny_global*nbits/8
+
+            ! use nbytes to compute RecSize.
+            ! this prevents integer overflow with large global grids using nbits
+            ! nx*ny*nbits > 2^31 -1 (i.e., global grid 9000x7054x64)
+            nbytes = nbits/bits_per_byte
+            RecSize = nx_global*ny_global*nbytes
+
             if (present(algn)) then
               ! If data is keept in blocks using given sizes (=algn)
               !  Used in eg. HYCOM binary files, which are stored as "blocks" dividable by 16384 bit (=algn)
@@ -131,6 +142,8 @@ subroutine ice_open_ext(nu, filename, nbits)
       integer (kind=int_kind), intent(in) :: &
            nu        , & ! unit number
            nbits         ! no. of bits per variable (0 for sequential access)
+
+      integer (kind=int_kind) :: RecSize, nbytes
 
       character (*) :: filename
 
@@ -150,7 +163,12 @@ subroutine ice_open_ext(nu, filename, nbits)
             nx = nx_global + 2*nghost
             ny = ny_global + 2*nghost
 
-            open(nu,file=filename,recl=nx*ny*nbits/8, &
+            ! use nbytes to compute RecSize.
+            ! this prevents integer overflow with large global grids using nbits
+            ! nx*ny*nbits > 2^31 -1 (i.e., global grid 9000x7054x64)
+            nbytes = nbits/bits_per_byte
+            RecSize = nx*ny*nbytes
+            open(nu,file=filename,recl=RecSize, &
                   form='unformatted',access='direct')
          endif                   ! nbits = 0
 

cicecore/shared/ice_init_column.F90

@@ -254,8 +254,6 @@ subroutine init_shortwave
 
       allocate(ztrcr_sw(nbtrcr_sw, ncat))
 
-      !!$OMP PARALLEL DO PRIVATE(iblk,i,j,n,ilo,ihi,jlo,jhi,this_block, &
-      !!$OMP                     l_print_point,debug,ipoint)
       do iblk=1,nblocks
 
          ! Initialize
@@ -387,12 +385,18 @@ subroutine init_shortwave
               enddo
             endif
 
+         enddo ! i
+         enddo ! j
+
       !-----------------------------------------------------------------
       ! Aggregate albedos 
+      ! Match loop order in coupling_prep for same order of operations
       !-----------------------------------------------------------------
 
-            do n = 1, ncat
-
+         do n = 1, ncat
+         do j = jlo, jhi
+         do i = ilo, ihi
+
                if (aicen(i,j,n,iblk) > puny) then
 
                   alvdf(i,j,iblk) = alvdf(i,j,iblk) &
@@ -422,7 +426,12 @@ subroutine init_shortwave
 
                endif ! aicen > puny
 
-            enddo  ! ncat
+         enddo ! i
+         enddo ! j
+         enddo ! ncat
+
+         do j = 1, ny_block
+         do i = 1, nx_block
 
       !----------------------------------------------------------------
       ! Store grid box mean albedos and fluxes before scaling by aice
@@ -432,14 +441,14 @@ subroutine init_shortwave
             alidf_ai  (i,j,iblk) = alidf  (i,j,iblk)
             alvdr_ai  (i,j,iblk) = alvdr  (i,j,iblk)
             alidr_ai  (i,j,iblk) = alidr  (i,j,iblk)
-            
+
             ! for history averaging
 !echmod?            cszn = c0
 !echmod            if (coszen(i,j,iblk) > puny) cszn = c1
 !echmod            do n = 1, nstreams
 !echmod               albcnt(i,j,iblk,n) = albcnt(i,j,iblk,n) + cszn
 !echmod            enddo
-            
+
       !----------------------------------------------------------------
       ! Save net shortwave for scaling factor in scale_factor
       !----------------------------------------------------------------

configuration/scripts/ice_in

@@ -118,6 +118,9 @@
     e_ratio         = 2.
     basalstress     = .false.
     k1              = 8.
+    k2              = 15.
+    alphab          = 20.
+    threshold_hw    = 30.
     coriolis        = 'latitude'
     kridge          = 1
     ktransport      = 1

configuration/scripts/options/set_nml.jra55

@@ -11,7 +11,7 @@ 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'
+atm_data_type   = 'JRA55_gx1'
 precip_units    = 'mks'
 ocn_data_dir    = 'default'
 bgc_data_dir    = 'default'

configuration/scripts/options/set_nml.jra55_gx1

@@ -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_gx1'
+precip_units    = 'mks'
+ocn_data_dir    = 'default'
+bgc_data_dir    = 'default'

configuration/scripts/options/set_nml.jra55_gx1_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_gx1'
+precip_units    = 'mks'
+ocn_data_dir    = 'default'
+bgc_data_dir    = 'default'

configuration/scripts/options/set_nml.jra55_gx3

@@ -0,0 +1,17 @@
+year_init      = 2005
+ice_ic  = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx3/iced_gx3_v5.nc'
+grid_file   = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/grid_gx3.bin'
+kmt_file    = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/kmt_gx3.bin'
+bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/global_gx3.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/gx3/JRA55'
+atm_data_format = 'nc'
+atm_data_type   = 'JRA55_gx3'
+precip_units    = 'mks'
+ocn_data_dir    = 'default'
+bgc_data_dir    = 'default'

configuration/scripts/options/set_nml.jra55_gx3_2008

@@ -0,0 +1,17 @@
+year_init      = 2008
+ice_ic  = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx3/iced_gx3_v5.nc'
+grid_file   = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/grid_gx3.bin'
+kmt_file    = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/kmt_gx3.bin'
+bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/global_gx3.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/gx3/JRA55'
+atm_data_format = 'nc'
+atm_data_type   = 'JRA55_gx3'
+precip_units    = 'mks'
+ocn_data_dir    = 'default'
+bgc_data_dir    = 'default'