Allow for read of tlat, tlon, anglet with popgrid (CICE-Consortium#463)

1/ To fix bug in ice_Transport driver .F90 (original pull request)

2/ Fix bug when cpp flag NEMO_IN_CICE is enabled. ice_step_mod.F90 uses variable raice without declaring.

3/ Allow for reading angle, tlon tlat if angle exist in netcdf grid file. If these are read they will not be calculated. This will cause differences when either of these three variables are used or dumped to output. If angle do not exist the routine will calculate angle, tlon and tlat as previous.

cicecore/cicedynB/dynamics/ice_transport_driver.F90

@@ -609,8 +609,8 @@ subroutine transport_remap (dt)
                                       asum_init(0), asum_final(0))
 
             if (l_stop) then
-               write (nu_diag,*) 'istep1, my_task, iblk =',     &
-                                  istep1, my_task, iblk
+               write (nu_diag,*) 'istep1, my_task =',     &
+                                  istep1, my_task
                write (nu_diag,*) 'transport: conservation error, cat 0'
                call abort_ice(subname//'ERROR: conservation error1')
             endif
@@ -623,8 +623,8 @@ subroutine transport_remap (dt)
                                       atsum_init(:,n), atsum_final(:,n))
 
                if (l_stop) then
-                  write (nu_diag,*) 'istep1, my_task, iblk, cat =',     &
-                                     istep1, my_task, iblk, n
+                  write (nu_diag,*) 'istep1, my_task, cat =',     &
+                                     istep1, my_task, n
                   write (nu_diag,*) 'transport: conservation error, cat ',n
                   call abort_ice(subname//'ERROR: conservation error2')
                endif

cicecore/cicedynB/general/ice_step_mod.F90

@@ -182,15 +182,17 @@ subroutine step_therm1 (dt, iblk)
       logical (kind=log_kind) :: & 
          prescribed_ice ! if .true., use prescribed ice instead of computed
 #endif
-
       real (kind=dbl_kind), intent(in) :: &
          dt      ! time step
 
       integer (kind=int_kind), intent(in) :: &
          iblk    ! block index
 
       ! local variables
-
+#ifdef CICE_IN_NEMO
+      real (kind=dbl_kind)    :: & 
+         raice              ! temporary reverse ice concentration
+#endif
       integer (kind=int_kind) :: &
          ilo,ihi,jlo,jhi, & ! beginning and end of physical domain
          i, j           , & ! horizontal indices

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -133,6 +133,10 @@ module ice_grid
       real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
          rndex_global       ! global index for local subdomain (dbl)
 
+      logical (kind=log_kind), private :: &
+         l_readCenter ! If anglet exist in grid file read it otherwise calculate it
+
+
 !=======================================================================
 
       contains
@@ -332,7 +336,6 @@ subroutine init_grid2
       real (kind=dbl_kind) :: &
          angle_0, angle_w, angle_s, angle_sw, &
          pi, pi2, puny
-!      real (kind=dbl_kind) :: ANGLET_dum
       logical (kind=log_kind), dimension(nx_block,ny_block,max_blocks):: &
          out_of_range
 
@@ -470,11 +473,10 @@ subroutine init_grid2
       !-----------------------------------------------------------------
       ! Compute ANGLE on T-grid
       !-----------------------------------------------------------------
-      ANGLET = c0
-
       if (trim(grid_type) == 'cpom_grid') then
          ANGLET(:,:,:) = ANGLE(:,:,:)
-      else
+      else if (.not. (l_readCenter)) then
+          ANGLET = c0
 
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
       !$OMP                     angle_0,angle_w,angle_s,angle_sw)
@@ -504,7 +506,8 @@ subroutine init_grid2
       enddo
       !$OMP END PARALLEL DO
       endif ! cpom_grid
-      if (trim(grid_type) == 'regional') then
+      if (trim(grid_type) == 'regional' .and. &
+          (.not. (l_readCenter))) then
          ! for W boundary extrapolate from interior
          !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
          do iblk = 1, nblocks
@@ -531,9 +534,9 @@ subroutine init_grid2
       call ice_timer_stop(timer_bound)
 
       call makemask          ! velocity mask, hemisphere masks
-
-      call Tlatlon           ! get lat, lon on the T grid
-
+      if (.not. (l_readCenter)) then
+         call Tlatlon           ! get lat, lon on the T grid
+      endif
       !-----------------------------------------------------------------
       ! bathymetry
       !-----------------------------------------------------------------
@@ -716,6 +719,7 @@ subroutine popgrid_nc
           field_loc_center, field_loc_NEcorner, &
           field_type_scalar, field_type_angle
       use ice_domain_size, only: max_blocks
+      use netcdf
 
       integer (kind=int_kind) :: &
          i, j, iblk, &
@@ -739,6 +743,12 @@ subroutine popgrid_nc
 
       type (block) :: &
          this_block           ! block information for current block
+
+      integer(kind=int_kind) :: &
+         varid
+      integer (kind=int_kind) :: &
+         status                ! status flag
+
 
       character(len=*), parameter :: subname = '(popgrid_nc)'
 
@@ -751,7 +761,7 @@ subroutine popgrid_nc
       call ice_open_nc(kmt_file,fid_kmt)
 
       diag = .true.       ! write diagnostic info
-
+      l_readCenter = .false.
       !-----------------------------------------------------------------
       ! topography
       !-----------------------------------------------------------------
@@ -806,11 +816,37 @@ subroutine popgrid_nc
       call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ANGLE
       call scatter_global(ANGLE, work_g1, master_task, distrb_info, &
                           field_loc_NEcorner, field_type_angle)
-
       ! fix ANGLE: roundoff error due to single precision
       where (ANGLE >  pi) ANGLE =  pi
       where (ANGLE < -pi) ANGLE = -pi
 
+      ! if grid file includes anglet then read instead
+      fieldname='anglet'
+      if (my_task == master_task) then
+         status = nf90_inq_varid(fid_grid, trim(fieldname) , varid)
+         if (status /= nf90_noerr) then
+            write(nu_diag,*) subname//' CICE will calculate angleT, TLON and TLAT'
+         else
+            write(nu_diag,*) subname//' angleT, TLON and TLAT is read from grid file'
+            l_readCenter = .true.
+         endif
+      endif
+      call broadcast_scalar(l_readCenter,master_task)
+      if (l_readCenter) then
+         call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) 
+         call scatter_global(ANGLET, work_g1, master_task, distrb_info, &
+                             field_loc_center, field_type_angle)
+         where (ANGLET >  pi) ANGLET =  pi
+         where (ANGLET < -pi) ANGLET = -pi
+         fieldname="tlon"
+         call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag)
+         call scatter_global(TLON, work_g1, master_task, distrb_info, &
+                             field_loc_center, field_type_scalar)
+         fieldname="tlat"
+         call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag)
+         call scatter_global(TLAT, work_g1, master_task, distrb_info, &
+                             field_loc_center, field_type_scalar)
+      endif
       !-----------------------------------------------------------------
       ! cell dimensions
       ! calculate derived quantities from global arrays to preserve 
@@ -820,7 +856,6 @@ subroutine popgrid_nc
       fieldname='htn'
       call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! HTN
       call primary_grid_lengths_HTN(work_g1)                  ! dxu, dxt
-
       fieldname='hte'
       call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! HTE
       call primary_grid_lengths_HTE(work_g1)                  ! dyu, dyt
@@ -831,7 +866,6 @@ subroutine popgrid_nc
          call ice_close_nc(fid_grid)
          call ice_close_nc(fid_kmt)
       endif
-
 #endif
       end subroutine popgrid_nc
 
@@ -1737,7 +1771,6 @@ subroutine Tlatlon
          enddo                  ! j         
       enddo                     ! iblk
       !$OMP END PARALLEL DO
-
       if (trim(grid_type) == 'regional') then
          ! for W boundary extrapolate from interior
          !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
@@ -2463,7 +2496,6 @@ subroutine read_basalstress_bathy
 
       ! use module
       use ice_read_write
-      use ice_communicate, only: my_task, master_task
       use ice_constants, only: field_loc_center, field_type_scalar
 
       ! local variables
@@ -2491,7 +2523,6 @@ subroutine read_basalstress_bathy
 
       if (my_task == master_task) then
          write(nu_diag,*) 'reading ',TRIM(fieldname)
-         write(*,*) 'reading ',TRIM(fieldname)
          call icepack_warnings_flush(nu_diag)
       endif
       call ice_read_nc(fid_init,1,fieldname,bathymetry,diag, &