Update CICE with CMIP changes.

cicecore/cicedynB/dynamics/ice_dyn_eap.F90

@@ -94,9 +94,9 @@ subroutine eap (dt)
           basal_stress_coeff, basalstress
       use ice_flux, only: rdg_conv, strairxT, strairyT, &
           strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, &
-          strtltx, strtlty, strocnx, strocny, strintx, strinty, &
-          strocnxT, strocnyT, &
-          Cbu, taubx, tauby, hwater, &
+          strtltx, strtlty, strocnx, strocny, strintx, strinty, taubx, tauby, &
+          strocnxT, strocnyT, strax, stray, &
+          Tbu, hwater, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
           stress12_1, stress12_2, stress12_3, stress12_4

cicecore/cicedynB/dynamics/ice_dyn_evp.F90

@@ -81,9 +81,9 @@ subroutine evp (dt)
       use ice_domain_size, only: max_blocks, ncat
       use ice_flux, only: rdg_conv, rdg_shear, strairxT, strairyT, &
           strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, &
-          strtltx, strtlty, strocnx, strocny, strintx, strinty, &
-          strocnxT, strocnyT, &
-          Cbu, taubx, tauby, hwater, &
+          strtltx, strtlty, strocnx, strocny, strintx, strinty, taubx, tauby, &
+          strocnxT, strocnyT, strax, stray, &
+          Tbu, hwater, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
           stress12_1, stress12_2, stress12_3, stress12_4

cicecore/cicedynB/general/ice_flux.F90

@@ -329,7 +329,6 @@ module ice_flux
       subroutine init_coupler_flux
 
       use ice_arrays_column, only: Cdn_atm
-      use ice_constants, only: p001
       use ice_flux_bgc, only: flux_bio_atm, flux_bio, faero_atm, &
            fnit, famm, fsil, fdmsp, fdms, fhum, fdust, falgalN, &
            fdoc, fdon, fdic, ffed, ffep
@@ -442,6 +441,7 @@ subroutine init_coupler_flux
       uocn  (:,:,:) = c0              ! surface ocean currents (m/s)
       vocn  (:,:,:) = c0
       frzmlt(:,:,:) = c0              ! freezing/melting potential (W/m^2)
+      frzmlt_init(:,:,:) = c0         ! freezing/melting potential (W/m^2)
       sss   (:,:,:) = 34.0_dbl_kind   ! sea surface salinity (ppt)
 
       do iblk = 1, size(Tf,3)
@@ -470,6 +470,7 @@ subroutine init_coupler_flux
       fsens   (:,:,:) = c0
       flat    (:,:,:) = c0
       fswabs  (:,:,:) = c0
+      fswint_ai(:,:,:) = c0
       flwout  (:,:,:) = -stefan_boltzmann*Tffresh**4   
                         ! in case atm model diagnoses Tsfc from flwout
       evap    (:,:,:) = c0
@@ -489,6 +490,7 @@ subroutine init_coupler_flux
       strocnyT(:,:,:) = c0    ! ice-ocean stress, y-direction (T-cell)
       fresh   (:,:,:) = c0
       fsalt   (:,:,:) = c0
+      fpond   (:,:,:) = c0
       fhocn   (:,:,:) = c0
       fswthru (:,:,:) = c0
       fresh_da(:,:,:) = c0    ! data assimilation
@@ -514,6 +516,7 @@ subroutine init_coupler_flux
 
       coszen  (:,:,:) = c0            ! Cosine of the zenith angle
       fsw     (:,:,:) = c0            ! shortwave radiation (W/m^2)
+      fswfac  (:,:,:) = c0
       scale_factor(:,:,:) = c1        ! shortwave scaling factor 
       wind    (:,:,:) = sqrt(uatm(:,:,:)**2 &
                            + vatm(:,:,:)**2)  ! wind speed, (m/s)
@@ -619,6 +622,7 @@ subroutine init_history_therm
       fsurf  (:,:,:) = c0
       fcondtop(:,:,:)= c0
       congel (:,:,:) = c0
+      fbot   (:,:,:) = c0
       frazil (:,:,:) = c0
       snoice (:,:,:) = c0
       dsnow  (:,:,:) = c0
@@ -637,7 +641,6 @@ subroutine init_history_therm
       fcondtopn (:,:,:,:) = c0
       flatn     (:,:,:,:) = c0
       fsensn    (:,:,:,:) = c0
-      fpond     (:,:,:) = c0
       fresh_ai  (:,:,:) = c0
       fsalt_ai  (:,:,:) = c0
       fhocn_ai  (:,:,:) = c0
@@ -654,10 +657,10 @@ subroutine init_history_therm
       Cdn_ocn(:,:,:) = dragio
       Cdn_atm(:,:,:) = (vonkar/log(zref/iceruf)) &
                      * (vonkar/log(zref/iceruf)) ! atmo drag for RASM
+      Cdn_atm_ratio(:,:,:)= c0
 
       if (formdrag) then
         Cdn_atm_rdg (:,:,:) = c0
-        Cdn_atm_ratio(:,:,:)= c0
         Cdn_atm_floe(:,:,:) = c0
         Cdn_atm_pond(:,:,:) = c0
         Cdn_atm_skin(:,:,:) = c0
@@ -685,7 +688,7 @@ end subroutine init_history_therm
 
       subroutine init_history_dyn
 
-      use ice_state, only: aice, vice, trcr
+      use ice_state, only: aice, vice, trcr, strength
 
       logical (kind=log_kind) :: &
           tr_iage
@@ -703,6 +706,7 @@ subroutine init_history_dyn
       sig2    (:,:,:) = c0
       taubx   (:,:,:) = c0
       tauby   (:,:,:) = c0
+      strength (:,:,:) = c0
       strocnx (:,:,:) = c0
       strocny (:,:,:) = c0
       strairx (:,:,:) = c0

cicecore/cicedynB/general/ice_forcing.F90

@@ -164,6 +164,21 @@ subroutine init_forcing_atmo
 
       character(len=*), parameter :: subname = '(init_forcing_atmo)'
 
+      fsw_data(:,:,:,:) = c0
+      cldf_data(:,:,:,:) = c0
+      fsnow_data(:,:,:,:) = c0
+      Tair_data(:,:,:,:) = c0
+      uatm_data(:,:,:,:) = c0
+      vatm_data(:,:,:,:) = c0
+      wind_data(:,:,:,:) = c0
+      strax_data(:,:,:,:) = c0
+      stray_data(:,:,:,:) = c0
+      Qa_data(:,:,:,:) = c0
+      rhoa_data(:,:,:,:) = c0
+      flw_data(:,:,:,:) = c0
+      sublim_data(:,:,:,:) = c0
+      frain_data(:,:,:,:) = c0
+
       fyear       = fyear_init + mod(nyr-1,ycycle) ! current year
       fyear_final = fyear_init + ycycle - 1 ! last year in forcing cycle
 
@@ -253,6 +268,11 @@ subroutine init_forcing_ocn(dt)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
 
+      sst_data(:,:,:,:) = c0
+      sss_data(:,:,:,:) = c0
+      uocn_data(:,:,:,:) = c0
+      vocn_data(:,:,:,:) = c0
+
       nbits = 64              ! double precision data
 
       if (restore_sst .or. restore_bgc) then
@@ -610,7 +630,7 @@ subroutine read_data (flag, recd, yr, ixm, ixx, ixp, &
          maxrec                  ! maximum record value
 
       real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks), &
-         intent(out) :: &
+         intent(inout) :: &
          field_data              ! 2 values needed for interpolation
 
       integer (kind=int_kind), intent(in) :: &
@@ -909,7 +929,7 @@ subroutine read_clim_data (readflag, recd, ixm, ixx, ixp, &
            field_type        ! type of field (scalar, vector, angle)
 
       real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks), &
-        intent(out) :: &
+        intent(inout) :: &
         field_data         ! 2 values needed for interpolation
 
       ! local variables
@@ -1257,7 +1277,6 @@ subroutine prepare_forcing (nx_block, ny_block, &
          ilo,ihi,jlo,jhi       ! beginning and end of physical domain
 
       real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
-         Tair    , & ! air temperature  (K)
          ANGLET  , & ! ANGLE converted to T-cells
          Tsfc    , & ! ice skin temperature
          sst     , & ! sea surface temperature
@@ -1270,6 +1289,7 @@ subroutine prepare_forcing (nx_block, ny_block, &
          cldf    , & ! cloud fraction
          frain   , & ! rainfall rate (kg/m^2 s)
          fsnow   , & ! snowfall rate (kg/m^2 s)
+         Tair    , & ! air temperature  (K)
          Qa      , & ! specific humidity (kg/kg)
          rhoa    , & ! air density (kg/m^3)
          uatm    , & ! wind velocity components (m/s)
@@ -1291,13 +1311,13 @@ subroutine prepare_forcing (nx_block, ny_block, &
          i, j
 
       real (kind=dbl_kind) :: workx, worky, &
-         precip_factor, zlvl0, secday, Tffresh
+         precip_factor, zlvl0, secday, Tffresh, puny
 
       logical (kind=log_kind) :: calc_strair
 
       character(len=*), parameter :: subname = '(prepare_forcing)'
 
-      call icepack_query_parameters(Tffresh_out=Tffresh)
+      call icepack_query_parameters(Tffresh_out=Tffresh, puny_out=puny)
       call icepack_query_parameters(secday_out=secday)
       call icepack_query_parameters(calc_strair_out=calc_strair)
       call icepack_warnings_flush(nu_diag)
@@ -1318,6 +1338,9 @@ subroutine prepare_forcing (nx_block, ny_block, &
          rhoa (i,j) = max(rhoa(i,j),c0)
          Qa   (i,j) = max(Qa(i,j),c0)
 
+         if (rhoa(i,j) .lt. puny) rhoa(i,j) = 1.3_dbl_kind            
+         if (Tair(i,j) .lt. puny) Tair(i,j) = Tffresh
+         if (Qa(i,j) .lt. puny) Qa(i,j) = 0.0035_dbl_kind
       enddo                     ! i
       enddo                     ! j
 
@@ -3635,6 +3658,8 @@ subroutine ocn_data_oned
       use ice_flux, only: sss, sst, Tf, uocn, vocn, ss_tltx, ss_tlty, &
             qdp, hmix, frzmlt
 
+      character(len=*), parameter :: subname = '(ocn_data_oned)'
+
       sss    (:,:,:) = 34.0_dbl_kind   ! sea surface salinity (ppt)
 
       call ocn_freezing_temperature
@@ -3796,7 +3821,7 @@ subroutine ocn_data_hadgem(dt)
      ! and change  units
      !----------------------------------------------------------------- 
 
-         !$OMP PARALLEL DO PRIVATE(iblk,i,j)
+         !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky)
          do iblk = 1, nblocks
             do j = 1, ny_block
             do i = 1, nx_block
@@ -3871,7 +3896,7 @@ subroutine read_data_nc_point (flag, recd, yr, ixm, ixx, ixp, &
            field_type        ! type of field (scalar, vector, angle)
 
       real (kind=dbl_kind), dimension(2), &
-         intent(out) :: &
+         intent(inout) :: &
          field_data              ! 2 values needed for interpolation
 
       character(len=*), parameter :: subname = '(read_data_nc_point)'
@@ -3886,6 +3911,8 @@ subroutine read_data_nc_point (flag, recd, yr, ixm, ixx, ixp, &
 
       call ice_timer_start(timer_readwrite)  ! reading/writing
 
+      field_data = c0 ! to satisfy intent(out) attribute
+
       if (istep1 > check_step) dbug = .true.  !! debugging
 
       if (my_task==master_task .and. (dbug)) then
@@ -4090,6 +4117,8 @@ subroutine ISPOL_data
 
       logical (kind=log_kind) :: read1
 
+      character(len=*), parameter :: subname = '(ISPOL_data)'
+
       call icepack_query_parameters(secday_out=secday)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &

cicecore/cicedynB/general/ice_step_mod.F90

@@ -75,6 +75,11 @@ subroutine prep_radiation (iblk)
 
       call ice_timer_start(timer_sw)      ! shortwave
 
+      alvdr_init(:,:,:) = c0
+      alvdf_init(:,:,:) = c0
+      alidr_init(:,:,:) = c0
+      alidf_init(:,:,:) = c0
+
          this_block = get_block(blocks_ice(iblk),iblk)         
          ilo = this_block%ilo
          ihi = this_block%ihi
@@ -142,7 +147,7 @@ subroutine step_therm1 (dt, iblk)
           fswthru, meltt, melts, meltb, congel, snoice, &
           flatn_f, fsensn_f, fsurfn_f, fcondtopn_f
       use ice_flux_bgc, only: dsnown, faero_atm, faero_ocn
-      use ice_grid, only: lmask_n, lmask_s
+      use ice_grid, only: lmask_n, lmask_s, tmask
       use ice_state, only: aice, aicen, aice_init, aicen_init, vicen_init, &
           vice, vicen, vsno, vsnon, trcrn, uvel, vvel, vsnon_init
 
@@ -267,6 +272,7 @@ subroutine step_therm1 (dt, iblk)
             enddo
          endif ! tr_aero
 
+         if (tmask(i,j,iblk)) &
          call icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero,                &
                             aicen_init  (i,j,:,iblk),                           &
                             vicen_init  (i,j,:,iblk), vsnon_init  (i,j,:,iblk), &
@@ -518,7 +524,7 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset)
       ! Aggregate the updated state variables (includes ghost cells). 
       !----------------------------------------------------------------- 
 
-         if (tmask(i,j,iblk)) &
+!        if (tmask(i,j,iblk)) &
          call icepack_aggregate (ncat,               aicen(i,j,:,iblk),   &
                                trcrn(i,j,1:ntrcr,:,iblk),               &
                                vicen(i,j,:,iblk), vsnon(i,j,  :,iblk),  &
@@ -906,6 +912,8 @@ subroutine ocean_mixed_layer (dt, iblk)
 
       use ice_arrays_column, only: Cdn_atm, Cdn_atm_ratio
       use ice_blocks, only: nx_block, ny_block
+      use ice_blocks, only: block, get_block
+      use ice_domain, only: blocks_ice
       use ice_flux, only: sst, Tf, Qa, uatm, vatm, wind, potT, rhoa, zlvl, &
            frzmlt, fhocn, fswthru, flw, flwout_ocn, fsens_ocn, flat_ocn, evap_ocn, &
            alvdr_ocn, alidr_ocn, alvdf_ocn, alidf_ocn, swidf, swvdf, swidr, swvdr, &
@@ -927,6 +935,7 @@ subroutine ocean_mixed_layer (dt, iblk)
          frzmlt_max = c1000   ! max magnitude of frzmlt (W/m^2)
 
       integer (kind=int_kind) :: &
+         ilo,ihi,jlo,jhi, & ! beginning and end of physical domain
          i, j           , & ! horizontal indices
          ij                 ! combined ij index
 
@@ -942,6 +951,8 @@ subroutine ocean_mixed_layer (dt, iblk)
       integer (kind=int_kind), dimension(nx_block*ny_block) :: &
          indxi, indxj    ! compressed indices for ocean cells
 
+      type (block) :: &
+         this_block         ! block information for current block
 
       !-----------------------------------------------------------------
 
@@ -958,8 +969,15 @@ subroutine ocean_mixed_layer (dt, iblk)
          icells = 0
          indxi(:) = 0
          indxj(:) = 0
-         do j = 1, ny_block
-         do i = 1, nx_block
+
+         this_block = get_block(blocks_ice(iblk),iblk)         
+         ilo = this_block%ilo
+         ihi = this_block%ihi
+         jlo = this_block%jlo
+         jhi = this_block%jhi
+
+         do j = jlo, jhi
+         do i = ilo, ihi
             if (tmask(i,j,iblk)) then
                icells = icells + 1
                indxi(icells) = i

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -302,6 +302,8 @@ subroutine init_grid2
       ! T-grid cell and U-grid cell quantities
       !-----------------------------------------------------------------
 
+      tarea(:,:,:) = c0
+
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
          this_block = get_block(blocks_ice(iblk),iblk)         
@@ -1489,6 +1491,8 @@ subroutine makemask
       !-----------------------------------------------------------------
 
       bm = c0
+      uvm = c0
+
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
          this_block = get_block(blocks_ice(iblk),iblk)         
@@ -2274,7 +2278,7 @@ subroutine read_basalstress_bathy
           write (nu_diag,*) ' '
           write (nu_diag,*) 'Initial ice file: ', trim(init_file)
           write (*,*) 'Initial ice file: ', trim(init_file)
-!         call flush(nu_diag)
+          call icepack_warnings_flush(nu_diag)
 
       endif
 
@@ -2285,7 +2289,7 @@ subroutine read_basalstress_bathy
       if (my_task == master_task) then
          write(nu_diag,*) 'reading ',TRIM(fieldname)
          write(*,*) 'reading ',TRIM(fieldname)
-!        call flush(nu_diag)
+         call icepack_warnings_flush(nu_diag)
       endif
       call ice_read_nc(fid_init,1,fieldname,bathymetry,diag, &
                     field_loc=field_loc_center, &
@@ -2295,7 +2299,7 @@ subroutine read_basalstress_bathy
 
       if (my_task == master_task) then
          write(nu_diag,*) 'closing file ',TRIM(init_file)
-!        call flush(nu_diag)
+         call icepack_warnings_flush(nu_diag)
       endif
 
       end subroutine read_basalstress_bathy