Update to CMIP7 diagnostics (#35)

* port CICE6 history features for CMIP variables: avg_ice_present, mask_ice_free_points
* add some error handling for variables not available in the ACCESS coupling
* add comment metadata for area and time weighting of variables
* use variable descriptions from data request, and set cell_methods
* add variables which are high priority, available and not previously included
* for snow variables, add an _intensive version, see https://airtable.com/appGoMZhCpkGESoVk/tblpo5L8maBIGlM1B/viwNNzrqK5oPL7zk2?blocks=hide

This change is based on the CMIP7 data request v1.2.2.2 - https://airtable.com/appGoMZhCpkGESoVk/tblpo5L8maBIGlM1B/viwNNzrqK5oPL7zk2?blocks=hide

Co-authored-by: Spencer Wong <88933912+blimlim@users.noreply.github.com>
Co-authored-by: SIobhan O'Farrell <spo599@gadi-login-05.gadi.nci.org.au>

Author: 3 people

drivers/access/CICE_RunMod.F90

@@ -62,14 +62,13 @@ subroutine CICE_Run
       use ice_timers, only: ice_timer_start, &
           ice_timer_stop, timer_couple, timer_step, &
           timer_from_atm, timer_into_atm, timer_from_ocn, timer_into_ocn
-      use ice_grid, only: t2ugrid_vector, u2tgrid_vector
+      use ice_grid, only: t2ugrid_vector
       integer (kind=int_kind) :: time_sec, itap, icpl_ai, tmp_time
       integer (kind=int_kind) :: rtimestamp_ai, stimestamp_ai
       integer (kind=int_kind) :: rtimestamp_io, stimestamp_io
       !receive and send timestamps (seconds)
       integer (kind=int_kind) :: imon 
 
-        logical ::  write_tmp_dump = .true.
 #endif
 
    !--------------------------------------------------------------------
@@ -83,9 +82,12 @@ subroutine CICE_Run
    !--------------------------------------------------------------------
 
 #ifdef ACCESS
-      write(il_out,*)'A <==> I coupling num_cpl_ai =          ',num_cpl_ai
-      write(il_out,*)' ice steps per ai interval num_ice_ai = ',num_ice_ai
-      write(il_out,*)' runtime, runtime0 =                    ',runtime, runtime0
+
+      if (my_task == master_task) then
+         write(il_out,*)'A <==> I coupling num_cpl_ai =          ',num_cpl_ai
+         write(il_out,*)' ice steps per ai interval num_ice_ai = ',num_ice_ai
+         write(il_out,*)' runtime, runtime0 =                    ',runtime, runtime0
+      endif
 
       time_sec = 0
 
@@ -94,7 +96,7 @@ subroutine CICE_Run
         !receive a2i fields 
         rtimestamp_ai = time_sec
         !call ice_timer_start(timer_from_atm)  ! atm/ice coupling
-        write(il_out,*)' calling from_atm at icpl_ai, time_sec = ', icpl_ai, time_sec
+      !   write(il_out,*)' calling from_atm at icpl_ai, time_sec = ', icpl_ai, time_sec
         !===========================
         call from_atm(rtimestamp_ai)
         !===========================
@@ -116,8 +118,8 @@ subroutine CICE_Run
           call t2ugrid_vector(io_strsu)
           call t2ugrid_vector(io_strsv)
 
-          write(il_out,'(a,3i10)') &
-                ' calling into_ocn at icpl_ai, itap, time_sec = ', icpl_ai, itap, time_sec
+         !  write(il_out,'(a,3i10)') &
+               !  ' calling into_ocn at icpl_ai, itap, time_sec = ', icpl_ai, itap, time_sec
           !call ice_timer_start(timer_into_ocn)  ! atm/ocn coupling
           !===========================
           !call check_iceberg_fields('chk_iceberg_i2o.nc')
@@ -161,8 +163,8 @@ subroutine CICE_Run
 
             stimestamp_ai = time_sec 
 
-            write(il_out,'(a,3i10)') &
-                ' calling into_atm at icpl_ai, itap, time_sec = ',icpl_ai, itap, time_sec
+            ! write(il_out,'(a,3i10)') &
+               !  ' calling into_atm at icpl_ai, itap, time_sec = ',icpl_ai, itap, time_sec
             !===========================
             call into_atm(stimestamp_ai)
             !===========================
@@ -292,6 +294,10 @@ subroutine ice_step
       use ice_algae, only: bgc_diags, write_restart_bgc
       use ice_zbgc, only: init_history_bgc, biogeochemistry
       use ice_zbgc_shared, only: skl_bgc
+#ifdef ACCESS
+      use ice_state, only: vsno, aice, tr_pond
+      use ice_flux, only: snowfrac
+#endif
 
       integer (kind=int_kind) :: &
          iblk        , & ! block index 
@@ -428,20 +434,12 @@ subroutine coupling_prep (iblk)
       use ice_constants, only: c0, c1, puny, rhofresh
       use ice_coupling, only: top_layer_Tandk_run, sfcflux_to_ocn
       use ice_domain_size, only: ncat
-      use ice_flux, only: alvdf, alidf, alvdr, alidr, albice, albsno, &
-          albpnd, albcnt, apeff_ai, coszen, fpond, fresh, &
-          alvdf_ai, alidf_ai, alvdr_ai, alidr_ai, fhocn_ai, &
-          fresh_ai, fsalt_ai, fsalt, &
-          fswthru_ai, fhocn, fswthru, scale_factor, &
-          swvdr, swidr, swvdf, swidf, Tf, Tair, Qa, strairxT, strairyt, &
-          fsens, flat, fswabs, flwout, evap, Tref, Qref, faero_ocn, &
-          fsurfn_f, flatn_f, scale_fluxes, frzmlt_init, frzmlt, &
-          snowfrac, snowfracn, evap_ice, evap_snow
+      use ice_flux
       use ice_grid, only: tmask
       use ice_ocean, only: oceanmixed_ice, ocean_mixed_layer
       use ice_shortwave, only: alvdfn, alidfn, alvdrn, alidrn, &
                                albicen, albsnon, albpndn, apeffn
-      use ice_state, only: aicen, aice, aice_init, nbtrcr
+      use ice_state, only: aicen, aice, aice_init, nbtrcr, tr_pond, vsno
       use ice_therm_shared, only: calc_Tsfc, heat_capacity
       use ice_timers, only: timer_couple, ice_timer_start, ice_timer_stop
       use ice_zbgc_shared, only: flux_bio, flux_bio_ai
@@ -527,8 +525,19 @@ subroutine coupling_prep (iblk)
 
             apeff_ai(i,j,iblk) = apeff_ai(i,j,iblk) &       ! for history
                + apeffn(i,j,n,iblk)*aicen(i,j,n,iblk)
-            snowfrac(i,j,iblk) = snowfrac(i,j,iblk) &       ! for history
-               + snowfracn(i,j,n,iblk)*aicen(i,j,n,iblk)
+
+            if ( .not. tr_pond .and. .not. calc_Tsfc ) then
+               ! calculate a snowfrac diagnostic in the same way the UM does
+               ! set snow fraction using JULES empirical formula based
+               ! on snow volume
+               ! ref: https://github.com/ACCESS-NRI/UM7/blob/6602dadd15c190ee37c6644190f52d428bc66917/umbase_hg3/src/atmosphere/short_wave_radiation/ftsa.F90#L201-L202
+               if (aice(i,j,iblk) > 2e-4) & 
+                  snowfrac(i,j,iblk) = c1 - exp(-p2*rhos*(vsno(i,j,iblk) / aice(i,j,iblk)))
+            else
+               snowfrac(i,j,iblk) = snowfrac(i,j,iblk) &       ! for history
+                  + snowfracn(i,j,n,iblk)*aicen(i,j,n,iblk)
+            endif
+
          enddo
          enddo
          enddo
@@ -555,6 +564,15 @@ subroutine coupling_prep (iblk)
             fsalt_ai  (i,j,iblk) = fsalt  (i,j,iblk)
             fhocn_ai  (i,j,iblk) = fhocn  (i,j,iblk)
             fswthru_ai(i,j,iblk) = fswthru(i,j,iblk)
+            fsens_ai  (i,j,iblk) = fsens(i,j,iblk)
+            flat_ai   (i,j,iblk) = flat(i,j,iblk)
+            fswabs_ai (i,j,iblk) = fswabs(i,j,iblk)
+            flwout_ai (i,j,iblk) = flwout(i,j,iblk)
+            evap_ai   (i,j,iblk) = evap(i,j,iblk) 
+            evap_ice_ai(i,j,iblk) = evap_ice(i,j,iblk) 
+            evap_snow_ai(i,j,iblk) = evap_snow(i,j,iblk) 
+            fcondtop_ai(i,j,iblk) = fcondtop(i,j,iblk) 
+            fsurf_ai(i,j,iblk) = fsurf(i,j,iblk) 
 
             if (nbtrcr > 0) then
             do k = 1, nbtrcr

drivers/access/cpl_forcing_handler.F90

@@ -424,9 +424,9 @@ subroutine get_lice_discharge(fname)
     case (4); myvar = 'FICEBERG_GC3_AVE'
   end select
   write(il_out,*)'(get_lice_discharge), iceberg = ', iceberg
-  write(il_out,'(a,a)') '(get_lice_discharge) reading in iceberg data, myvar= ',trim(myvar)
+  ! write(il_out,'(a,a)') '(get_lice_discharge) reading in iceberg data, myvar= ',trim(myvar)
   do im = 1, 12
-    write(il_out,*) '(get_lice_discharge) reading in data, month= ',im
+    ! write(il_out,*) '(get_lice_discharge) reading in data, month= ',im
     call ice_read_nc(ncid_i2o, im, trim(myvar), vwork, dbug)
 
     ! Restrict iceberg fluxes to ocean points

drivers/access/cpl_interface.F90

@@ -92,24 +92,23 @@ subroutine prism_init
   ! Initialise MPI
   mpiflag = .FALSE.
   call MPI_Initialized (mpiflag, ierror)
-  print *,'CICE: (prism_init) BF MPI_INIT, mpiflag = ',mpiflag
+  ! print *,'CICE: (prism_init) BF MPI_INIT, mpiflag = ',mpiflag
 
   if ( .not. mpiflag ) then
     call MPI_INIT(ierror)
   endif
 
   call MPI_Initialized (mpiflag, ierror)
-  print *, 'CICE: (prism_init) AF MPI_INIT, mpiflag = ',mpiflag
+  ! print *, 'CICE: (prism_init) AF MPI_INIT, mpiflag = ',mpiflag
 
-  print *
-  print *, 'CICE: (prism_init) calling prism_init_comp_proto...'
+  ! print *, 'CICE: (prism_init) calling prism_init_comp_proto...'
 
   call prism_init_comp_proto (il_comp_id, cp_modnam, ierror)
 
   if (ierror /= PRISM_Ok) then 
       call prism_abort_proto(il_comp_id, 'cice prism_init','STOP 1')
-  else
-      print *, 'CICE: (prism_init) called prism_init_comp_proto !'
+  ! else
+  !     if ( my_task == master_task ) print *, 'CICE: (prism_init) called prism_init_comp_proto !'
   endif
 
   !B: the following part may not be really needed(?)
@@ -132,8 +131,8 @@ subroutine prism_init
   if (ierror /= PRISM_Ok) then
       print *, 'CICE: (prism_init) Error in MPI_Buffer_Attach.'
       call prism_abort_proto(il_comp_id, 'cice prism_init','STOP 2')
-  else
-      print *, 'CICE: (prism_init) MPI_Buffer_Attach ok!'
+  ! else
+  !     print *, 'CICE: (prism_init) MPI_Buffer_Attach ok!'
   endif
   !
   ! PSMILe attribution of local communicator.
@@ -147,21 +146,21 @@ subroutine prism_init
       print *, 'CICE: Error in prism_get_localcomm_proto'
       call prism_abort_proto(il_comp_id, 'cice prism_init','STOP 3')
   else
-      print *, 'CICE: _get_localcomm_ OK! il_commlocal= ',il_commlocal
+      if ( my_task == master_task ) print *, 'CICE: _get_localcomm_ OK! il_commlocal= ',il_commlocal
   endif
 
   !
   ! Inquire if model is parallel or not and open the process log file 
   !
  ! print *, '* CICE: Entering init_cpl.....'
 
-  print *, '* CICE (prism_init) calling MPI_Comm_Size ...'
+  !print *, '* CICE (prism_init) calling MPI_Comm_Size ...'
   call MPI_Comm_Size(il_commlocal, il_nbtotproc, ierror)
-  print *, '* CICE (prism_init) calling MPI_Comm_Rank ...'
+  !print *, '* CICE (prism_init) calling MPI_Comm_Rank ...'
   call MPI_Comm_Rank(il_commlocal, my_task, ierror)
 
-  print *, '* CICE (prism_init) il_commlocal, il_nbtotproc, my_task = '
-  print *, '* CICE (prism_init) ', il_commlocal, il_nbtotproc, my_task
+  if ( my_task == master_task ) print *, '* CICE (prism_init) il_commlocal, il_nbtotproc, my_task = '
+  if ( my_task == master_task ) print *, '* CICE (prism_init) ', il_commlocal, il_nbtotproc, my_task
   !
   il_nbcplproc = il_nbtotproc   !multi-process coupling (real parallel cpl)!
   !il_nbcplproc = 1               !mono process coupling
@@ -258,18 +257,18 @@ subroutine init_cpl
 !     end do
 
 !!debug
-  if (my_task == 0) then
-  write(il_out,*) "all block info:"
-  do iblk=1,nblocks_tot
-      this_block = get_block(iblk,iblk)
-      ilo = this_block%ilo
-      ihi = this_block%ihi
-      jlo = this_block%jlo
-      jhi = this_block%jhi
-      write(il_out,*) '   this block: cpu, iblock, jblock=', distrb_info%blockLocation(iblk)-1, this_block%iblock, this_block%jblock
-      write(il_out,*) '   block:', iblk, "gilo, gjlo, gihi, gjhi=", this_block%i_glob(ilo), this_block%j_glob(jlo), this_block%i_glob(ihi), this_block%j_glob(jhi)
-  end do
-  end if
+  ! if (my_task == 0) then
+  ! write(il_out,*) "all block info:"
+  ! do iblk=1,nblocks_tot
+  !     this_block = get_block(iblk,iblk)
+  !     ilo = this_block%ilo
+  !     ihi = this_block%ihi
+  !     jlo = this_block%jlo
+  !     jhi = this_block%jhi
+  !     write(il_out,*) '   this block: cpu, iblock, jblock=', distrb_info%blockLocation(iblk)-1, this_block%iblock, this_block%jblock
+  !     write(il_out,*) '   block:', iblk, "gilo, gjlo, gihi, gjhi=", this_block%i_glob(ilo), this_block%j_glob(jlo), this_block%i_glob(ihi), this_block%j_glob(jhi)
+  ! end do
+  ! end if
 
   do iblk=1,nblocks_tot
 
@@ -781,8 +780,8 @@ subroutine from_atm(isteps)
       call ncheck( nf_open('fields_a2i_in_ice.nc',nf_write,ncid) )
       call write_nc_1Dtime(real(isteps),currstep,'time',ncid)
     endif
-    write(il_out,*)
-    write(il_out,*) '(from_atm) Total number of fields to be rcvd: ', nrecv_a2i
+    ! write(il_out,*)
+    ! write(il_out,*) '(from_atm) Total number of fields to be rcvd: ', nrecv_a2i
   endif
 
   if (debug) write(il_out,*) "prism_get from_atm at sec: ", isteps

io_netcdf/ice_history_write.F90

@@ -668,6 +668,11 @@ subroutine ice_hist_create(ns, ncfile, ncid, var, coord_var, var_nverts, var_nz)
                 call check(nf90_put_att(ncid,varid,'cell_measures', &
                             avail_hist_fields(n)%vcellmeas), &
                            'put att cell_measures '//avail_hist_fields(n)%vname)
+                if (avail_hist_fields(n)%vcomment /= "none") then
+                  call check(nf90_put_att(ncid,varid,'comment', &
+                            avail_hist_fields(n)%vcomment), &
+                            'put att comment '//avail_hist_fields(n)%vname)
+                endif
                 call check(nf90_put_att(ncid,varid,'missing_value',spval), &
                            'put att missing_value '//avail_hist_fields(n)%vname)
                 call check(nf90_put_att(ncid,varid,'_FillValue',spval), &
@@ -679,9 +684,20 @@ subroutine ice_hist_create(ns, ncfile, ncid, var, coord_var, var_nverts, var_nz)
             if (hist_avg) then
                 if (TRIM(avail_hist_fields(n)%vname)/='sig1' .or. &
                     TRIM(avail_hist_fields(n)%vname)/='sig2') then
-
-                    call check(nf90_put_att(ncid,varid,'cell_methods','time: mean'), &
+                    if (avail_hist_fields(n)%avg_ice_present) then
+                      call check(nf90_put_att(ncid,varid,'cell_methods',&
+                        'area: time: mean where sea_ice (mask=siconc)'), &
+                        'put att cell methods time mean '//avail_hist_fields(n)%vname)
+                    else
+                      if (TRIM(avail_hist_fields(n)%vname(1:2))/='si') then !native diags
+                        call check(nf90_put_att(ncid,varid,'cell_methods','time: mean'), &
+                          'put att cell methods time mean '//avail_hist_fields(n)%vname)
+                      else !cmip diags
+                        call check(nf90_put_att(ncid,varid,'cell_methods', &
+                          'area: mean where sea time: mean'), &
                           'put att cell methods time mean '//avail_hist_fields(n)%vname)
+                      endif
+                    endif
                 endif
             endif
 
@@ -739,6 +755,11 @@ subroutine ice_hist_create(ns, ncfile, ncid, var, coord_var, var_nverts, var_nz)
                         avail_hist_fields(n)%vcellmeas)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining cell measures for '//avail_hist_fields(n)%vname)
+            if (avail_hist_fields(n)%vcomment /= "none") then
+              call check(nf90_put_att(ncid,varid,'comment', &
+                        avail_hist_fields(n)%vcomment), &
+                        'put att comment '//avail_hist_fields(n)%vname)
+            endif
             status = nf90_put_att(ncid,varid,'missing_value',spval)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining missing_value for '//avail_hist_fields(n)%vname)
@@ -802,6 +823,11 @@ subroutine ice_hist_create(ns, ncfile, ncid, var, coord_var, var_nverts, var_nz)
                         avail_hist_fields(n)%vcellmeas)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining cell measures for '//avail_hist_fields(n)%vname)
+            if (avail_hist_fields(n)%vcomment /= "none") then
+              call check(nf90_put_att(ncid,varid,'comment', &
+                        avail_hist_fields(n)%vcomment), &
+                        'put att comment '//avail_hist_fields(n)%vname)
+            endif
             status = nf90_put_att(ncid,varid,'missing_value',spval)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining missing_value for '//avail_hist_fields(n)%vname)
@@ -851,6 +877,11 @@ subroutine ice_hist_create(ns, ncfile, ncid, var, coord_var, var_nverts, var_nz)
                         avail_hist_fields(n)%vcellmeas)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining cell measures for '//avail_hist_fields(n)%vname)
+            if (avail_hist_fields(n)%vcomment /= "none") then
+              call check(nf90_put_att(ncid,varid,'comment', &
+                        avail_hist_fields(n)%vcomment), &
+                        'put att comment '//avail_hist_fields(n)%vname)
+            endif
             status = nf90_put_att(ncid,varid,'missing_value',spval)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining missing_value for '//avail_hist_fields(n)%vname)
@@ -901,6 +932,11 @@ subroutine ice_hist_create(ns, ncfile, ncid, var, coord_var, var_nverts, var_nz)
                         avail_hist_fields(n)%vcellmeas)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining cell measures for '//avail_hist_fields(n)%vname)
+            if (avail_hist_fields(n)%vcomment /= "none") then
+              call check(nf90_put_att(ncid,varid,'comment', &
+                        avail_hist_fields(n)%vcomment), &
+                        'put att comment '//avail_hist_fields(n)%vname)
+            endif
             status = nf90_put_att(ncid,varid,'missing_value',spval)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining missing_value for '//avail_hist_fields(n)%vname)
@@ -966,6 +1002,11 @@ subroutine ice_hist_create(ns, ncfile, ncid, var, coord_var, var_nverts, var_nz)
                         avail_hist_fields(n)%vcellmeas)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining cell measures for '//avail_hist_fields(n)%vname)
+            if (avail_hist_fields(n)%vcomment /= "none") then
+              call check(nf90_put_att(ncid,varid,'comment', &
+                        avail_hist_fields(n)%vcomment), &
+                        'put att comment '//avail_hist_fields(n)%vname)
+            endif
             status = nf90_put_att(ncid,varid,'missing_value',spval)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining missing_value for '//avail_hist_fields(n)%vname)
@@ -1031,6 +1072,11 @@ subroutine ice_hist_create(ns, ncfile, ncid, var, coord_var, var_nverts, var_nz)
                         avail_hist_fields(n)%vcellmeas)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining cell measures for '//avail_hist_fields(n)%vname)
+            if (avail_hist_fields(n)%vcomment /= "none") then
+              call check(nf90_put_att(ncid,varid,'comment', &
+                        avail_hist_fields(n)%vcomment), &
+                        'put att comment '//avail_hist_fields(n)%vname)
+            endif
             status = nf90_put_att(ncid,varid,'missing_value',spval)
             if (status /= nf90_noerr) call abort_ice( &
                'Error defining missing_value for '//avail_hist_fields(n)%vname)
@@ -1094,9 +1140,10 @@ subroutine ice_hist_create(ns, ncfile, ncid, var, coord_var, var_nverts, var_nz)
         call check(nf90_put_att(ncid,nf90_global,'comment2',title), &
                    'global attribute comment2')
 
-        title = 'CF-1.0'
-        call check(nf90_put_att(ncid,nf90_global,'conventions',title), &
-                   'global attribute conventions')
+        ! TO-DO: Update output for CF compliance !
+        ! title = 'CF-1.0'
+        ! call check(nf90_put_att(ncid,nf90_global,'conventions',title), &
+        !            'global attribute conventions')
 
         call date_and_time(date=current_date, time=current_time)
         write(start_time,1000) current_date(1:4), current_date(5:6), &