+(*)Fix rotation of coupler_type variables

  This commit adds the necessary information about the turns of the model grid
relative to the (unrotated) coupler_type data fields that are inside of the
forcing type and surface_state type and are used with passive tracers, so that
certain tracer packages can now be used with rotated grids.  The previous
version had problems where the model would not run when ROTATE_INTEX = True and
the CFCs (and probably other passive tracers) were being used, as noted at
NOAA-GFDL/issues/621.  These problems have now been fixed.

  There are new calls to coupler_type_spawn() in allocate_forcing_by_ref() and
allocate_surface_state() that manage the creation of the coupler_2d_bt_types for
unrotated types based on information from their rotated counterparts.

  There is a new optional turns arguments to allocate_forcing_by_ref() and new
optional sfc_state_in and turns arguments to allocate_surface_state(), and these
are now being used in at least 6 places where unrotated temporary surface_state
or forcing types are being set up.

  There are also new optional turns argument to extract_coupler_type_data() and
set_coupler_type_data() that are used to deal with the fact that the internal
data arrays in the coupler_bc_types are never rotated, unlike the other MOM6
arrays, because they have to be passed directly into the generic tracer code.
These new turns arguments are used in 14 calls from various tracer packages,
including in 6 calls from the OCMIP2_CFC code.

  There are 4 new calls to deallocate_surface_state() or
deallocate_forcing_type() that were added to avoid (presumably minor) memory
leaks when grid rotation is enabled.  These new calls are in
initialize_ice_shelf_fluxes(), shelf_calc_flux() and in the surface flux
diagnostics section of step_MOM().

  All answers are bitwise identical in any cases that worked previously, but
some cases with grid rotation that previously were failing with cryptic error
messages are now running successfully.  There are new optional arguments to
several publicly visible routines

src/core/MOM.F90

@@ -624,7 +624,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
     call rotate_mech_forcing(forces_in, turns, forces)
 
     allocate(fluxes)
-    call allocate_forcing_type(fluxes_in, G, fluxes)
+    call allocate_forcing_type(fluxes_in, G, fluxes, turns=turns)
     call rotate_forcing(fluxes_in, fluxes, turns)
   else
     forces => forces_in
@@ -1044,6 +1044,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
 
   ! Do diagnostics that only occur at the end of a complete forcing step.
   if (cycle_end) then
+    if (showCallTree) call callTree_waypoint("Do cycle end diagnostics (step_MOM)")
     if (CS%rotate_index) then
       allocate(sfc_state_diag)
       call rotate_surface_state(sfc_state, sfc_state_diag, G, turns)
@@ -1063,6 +1064,10 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
     endif
     call disable_averaging(CS%diag)
     call cpu_clock_end(id_clock_diagnostics)
+    if (CS%rotate_index) then
+      call deallocate_surface_state(sfc_state_diag)
+    endif
+    if (showCallTree) call callTree_waypoint("Done with end cycle diagnostics (step_MOM)")
   endif
 
   ! Accumulate the surface fluxes for assessing conservation
@@ -3688,8 +3693,8 @@ subroutine extract_surface_state(CS, sfc_state_in)
   if (CS%rotate_index) then
     allocate(sfc_state)
     call allocate_surface_state(sfc_state, G, use_temperature, &
-         do_integrals=.true., omit_frazil=.not.associated(CS%tv%frazil),&
-         use_iceshelves=use_iceshelves)
+              do_integrals=.true., omit_frazil=.not.associated(CS%tv%frazil),&
+              use_iceshelves=use_iceshelves, sfc_state_in=sfc_state_in, turns=turns)
   else
     sfc_state => sfc_state_in
   endif

src/core/MOM_forcing_type.F90

@@ -6,7 +6,7 @@ module MOM_forcing_type
 use MOM_array_transform, only : rotate_array, rotate_vector, rotate_array_pair
 use MOM_coupler_types, only : coupler_2d_bc_type, coupler_type_destructor
 use MOM_coupler_types, only : coupler_type_increment_data, coupler_type_initialized
-use MOM_coupler_types, only : coupler_type_copy_data
+use MOM_coupler_types, only : coupler_type_copy_data, coupler_type_spawn
 use MOM_cpu_clock,     only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
 use MOM_debugging,     only : hchksum, uvchksum
 use MOM_diag_mediator, only : post_data, register_diag_field, register_scalar_field
@@ -2628,7 +2628,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
   if (turns /= 0) then
     G => diag%G
     allocate(fluxes)
-    call allocate_forcing_type(fluxes_in, G, fluxes)
+    call allocate_forcing_type(fluxes_in, G, fluxes, turns=turns)
     call rotate_forcing(fluxes_in, fluxes, turns)
   else
     G => G_in
@@ -3309,13 +3309,16 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
 end subroutine allocate_forcing_by_group
 
 !> Allocate elements of a new forcing type based on their status in an existing type.
-subroutine allocate_forcing_by_ref(fluxes_ref, G, fluxes)
+subroutine allocate_forcing_by_ref(fluxes_ref, G, fluxes, turns)
   type(forcing),         intent(in)  :: fluxes_ref !< Reference fluxes
   type(ocean_grid_type), intent(in)  :: G          !< Grid metric of target fluxes
   type(forcing),         intent(out) :: fluxes     !< Target fluxes
+  integer,     optional, intent(in)  :: turns      !< If present, the number of counterclockwise
+                                                   !! quarter turns to use on the new grid.
 
   logical :: do_ustar, do_taumag, do_water, do_heat, do_salt, do_press, do_shelf
   logical :: do_iceberg, do_heat_added, do_buoy
+  logical :: even_turns  ! True if turns is absent or even
 
   call get_forcing_groups(fluxes_ref, do_water, do_heat, do_ustar, do_taumag, do_press, &
       do_shelf, do_iceberg, do_salt, do_heat_added, do_buoy)
@@ -3354,6 +3357,19 @@ subroutine allocate_forcing_by_ref(fluxes_ref, G, fluxes)
   ! This flag would normally be set by a control flag in allocate_forcing_type.
   ! Here we copy the flag from the reference forcing.
   fluxes%gustless_accum_bug = fluxes_ref%gustless_accum_bug
+
+  if (coupler_type_initialized(fluxes_ref%tr_fluxes)) then
+    ! The data fields in the coupler_2d_bc_type are never rotated.
+    even_turns = .true. ; if (present(turns)) even_turns = (modulo(turns, 2) == 0)
+    if (even_turns) then
+      call coupler_type_spawn(fluxes_ref%tr_fluxes, fluxes%tr_fluxes, &
+                (/G%isc,G%isc,G%iec,G%iec/), (/G%jsc,G%jsc,G%jec,G%jec/))
+    else
+      call coupler_type_spawn(fluxes_ref%tr_fluxes, fluxes%tr_fluxes, &
+                (/G%jsc,G%jsc,G%jec,G%jec/), (/G%isc,G%isc,G%iec,G%iec/))
+    endif
+  endif
+
 end subroutine allocate_forcing_by_ref
 
 

src/core/MOM_variables.F90

@@ -338,14 +338,14 @@ module MOM_variables
 !! the ocean model. Unused fields are unallocated.
 subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
                                   gas_fields_ocn, use_meltpot, use_iceshelves, &
-                                  omit_frazil)
+                                  omit_frazil, sfc_state_in, turns)
   type(ocean_grid_type), intent(in)    :: G                !< ocean grid structure
   type(surface),         intent(inout) :: sfc_state        !< ocean surface state type to be allocated.
   logical,     optional, intent(in)    :: use_temperature  !< If true, allocate the space for thermodynamic variables.
   logical,     optional, intent(in)    :: do_integrals     !< If true, allocate the space for vertically
                                                            !! integrated fields.
   type(coupler_1d_bc_type), &
-               optional, intent(in)    :: gas_fields_ocn  !< If present, this type describes the ocean
+               optional, intent(in)    :: gas_fields_ocn   !< If present, this type describes the
                                               !! ocean and surface-ice fields that will participate
                                               !! in the calculation of additional gas or other
                                               !! tracer fluxes, and can be used to spawn related
@@ -355,9 +355,20 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
                                                            !! under ice shelves.
   logical,     optional, intent(in)    :: omit_frazil      !< If present and false, do not allocate the space to
                                                            !! pass frazil fluxes to the coupler
+  type(surface), &
+               optional, intent(in)    :: sfc_state_in     !< If present and its tr_fields are initialized,
+                                              !! this type describes the ocean and surface-ice fields that
+                                              !! will participate in the calculation of additional gas or
+                                              !! other tracer fluxes, and can be used to spawn related
+                                              !! internal variables in the ice model.  If gas_fields_ocn
+                                              !! is present, it is used and tr_fields_in is ignored.
+  integer,     optional, intent(in)    :: turns  !< If present, the number of counterclockwise quarter
+                                                 !! turns to use on the new grid.
 
   ! local variables
   logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil
+  logical :: even_turns  ! True if turns is absent or even
+  integer :: tr_field_i_mem(4), tr_field_j_mem(4)
   integer :: is, ie, js, je, isd, ied, jsd, jed
   integer :: isdB, iedB, jsdB, jedB
 
@@ -405,9 +416,22 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
     allocate(sfc_state%tauy_shelf(isd:ied,JsdB:JedB), source=0.0)
   endif
 
-  if (present(gas_fields_ocn)) &
+  ! The data fields in the coupler_2d_bc_type are never rotated.
+  even_turns = .true. ; if (present(turns)) even_turns = (modulo(turns, 2) == 0)
+  if (even_turns) then
+    tr_field_i_mem(1:4) = (/is,is,ie,ie/) ; tr_field_j_mem(1:4) = (/js,js,je,je/)
+  else
+    tr_field_i_mem(1:4) = (/js,js,je,je/) ; tr_field_j_mem(1:4) = (/is,is,ie,ie/)
+  endif
+  if (present(gas_fields_ocn)) then
     call coupler_type_spawn(gas_fields_ocn, sfc_state%tr_fields, &
-                            (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.)
+                            tr_field_i_mem, tr_field_j_mem, as_needed=.true.)
+  elseif (present(sfc_state_in)) then
+    if (coupler_type_initialized(sfc_state_in%tr_fields)) then
+      call coupler_type_spawn(sfc_state_in%tr_fields, sfc_state%tr_fields, &
+                              tr_field_i_mem, tr_field_j_mem, as_needed=.true.)
+    endif
+  endif
 
   sfc_state%arrays_allocated = .true.
 
@@ -438,10 +462,10 @@ end subroutine deallocate_surface_state
 
 !> Rotate the surface state fields from the input to the model indices.
 subroutine rotate_surface_state(sfc_state_in, sfc_state, G, turns)
-  type(surface), intent(in) :: sfc_state_in
-  type(surface), intent(inout) :: sfc_state
-  type(ocean_grid_type), intent(in) :: G
-  integer, intent(in) :: turns
+  type(surface), intent(in) :: sfc_state_in  !< The input unrotated surface state type that is the data source.
+  type(surface), intent(inout) :: sfc_state  !< The rotated surface state type whose arrays will be filled in
+  type(ocean_grid_type), intent(in) :: G     !< The ocean grid structure
+  integer, intent(in) :: turns   !< The number of counterclockwise quarter turns to use on the rotated grid.
 
   logical :: use_temperature, do_integrals, use_melt_potential, use_iceshelves
 
@@ -454,13 +478,9 @@ subroutine rotate_surface_state(sfc_state_in, sfc_state, G, turns)
       .and. allocated(sfc_state_in%tauy_shelf)
 
   if (.not. sfc_state%arrays_allocated) then
-    call allocate_surface_state(sfc_state, G, &
-        use_temperature=use_temperature, &
-        do_integrals=do_integrals, &
-        use_meltpot=use_melt_potential, &
-        use_iceshelves=use_iceshelves &
-    )
-    sfc_state%arrays_allocated = .true.
+    call allocate_surface_state(sfc_state, G, use_temperature=use_temperature, &
+            do_integrals=do_integrals, use_meltpot=use_melt_potential, &
+            use_iceshelves=use_iceshelves, sfc_state_in=sfc_state_in, turns=turns)
   endif
 
   if (use_temperature) then