IMPORTANT: Active and Passive tracer and Limiting bottom stress

[arango]

The following updates to the code were implemented:

• Added the switch LtracerCLM to manage the tracer climatology field to process for active and passive tracers in all the standard input files (ocean_*.in, biology, and sediment). For example, in ocean_*.in we have now:

  ! Logical switches (TRUE/FALSE) to specify which variables to process for
  ! tracers climatology: [1:NAT+NPT,Ngrids]. See glossary below for details.
  
    LtracerCLM == T T                        ! temperature, salinity, inert
  

  This will allow us to only process and operate on the desired tracer climatology field. In biological tracers, it will make sense to operate on the nutrient climatology and not on the higher trophic level model constituents, for example. The tracer climatology tclm array will only contains the desired climatology fields to save memory usage. They are processed as follows to read and interpolate:

        ic=0
        DO i=1,NT(ng)
          IF (LtracerCLM(i,ng)) THEN
            ic=ic+1
            CALL get_3dfld (ng, iNLM, idTclm(i), CLM(ng)%ncid,            &
       &                    1, CLM(ng), update(1),                        &
       &                    LBi, UBi, LBj, UBj, 1, N(ng), 2, 1,           &
  #   ifdef MASKING
       &                    GRID(ng) % rmask,                             &
  #   endif
       &                    CLIMA(ng) % tclmG(:,:,:,:,ic))
            IF (exit_flag.ne.NoError) RETURN
          END IF
        END DO
  
  ...
  
        ic=0
        DO itrc=1,NT(ng)
          IF (LtracerCLM(itrc,ng)) THEN
            ic=ic+1
            CALL set_3dfld_tile (ng, tile, iNLM, idTclm(itrc),            &
       &                         LBi, UBi, LBj, UBj, 1, N(ng),            &
       &                         CLIMA(ng)%tclmG(:,:,:,:,ic),             &
       &                         CLIMA(ng)%tclm (:,:,:,ic),               &
       &                         update)
            IF (exit_flag.ne.NoError) RETURN
          END IF
        END DO
  

  Notice that we have now the full tracer DO-loop.

In step3d_t.F, we now have:

# if defined TCLM_NUDGING && defined TCLIMATOLOGY
!
!  Nudge towards tracer climatology.
!
        IF (LtracerCLM(itrc,ng)) THEN
          ic=ic+1
          DO k=1,N(ng)
            DO j=JstrR,JendR
              DO i=IstrR,IendR
                t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+                  &
     &                             dt(ng)*Tnudgcof(i,j,itrc)*           &
     &                             (tclm(i,j,k,ic)-t(i,j,k,nnew,itrc))
              END DO
            END DO
          END DO
        END IF
# endif

• Added processing of all passive tracer quadratic terms: <t*t>, <u*t>, <v*t>, <Huon*t>, and <Hvom*t>. This will facilitate computing budgets for any passive tracer. The desired fields are processed and written into the average files according to the logical switches:

  Aout(idTTav) == F F     ! temp2, ...         quadratic <t*t> tracer terms
  Aout(idUTav) == F F     ! utemp, ...         quadratic <u*t> tracer terms
  Aout(idVTav) == F F     ! vtemp, ...         quadratic <v*t> tracer terms
  Aout(iHUTav) == F F     ! Huontemp, ...      tracer volume flux, <Huon*t>
  Aout(iHVTav) == F F     ! Hvomtemp, ...      tracer volume flux, <Hvom*t>
  

  This also required changes to all the standard input files (ocean_*.in, biology, and sediment).

Warming: The file varinfo.dat was also modified for the generic quadratic tracers metadata:

'tracer2'                                          ! Output
  'time-averaged squared tracer'
  'tracer2'                                        ! [tracer2]
  'tracer2, scalar, series'
  'ocean_time'
  'idTTav(itrc)'
  'r3dvar'
  1.0d0

'Huontracer'                                       ! Output
  'time-averaged tracer u-volume flux'
  'meter3 second-1 tracer'                         ! [m3/s tracer]
  'Huontracer, scalar, series'
  'ocean_time'
  'iHUTav(itrc)'
  'u3dvar'
 1.0d0

'utracer'                                          ! Output
  'time-averaged u-momentum times tracer'
  'meter second-1 tracer'                          ! [m/s tracer]
  'utracer, scalar, series'
  'ocean_time'
  'idUTav(itrc)'
  'u3dvar'
  1.0d0

'Hvomtracer'                                       ! Output
  'time-averaged tracer v-volume flux'
  'meter3 second-1 tracer'                         ! [m3/s tracer]
  'Hvomtracer, scalar, series'
  'ocean_time'
  'iHVTav(itrc)'
  'v3dvar'
 1.0d0

'vtracer'                                          ! Output
  'time-averaged v-momentum times tracer'
  'meter second-1 tracer'                          ! [m/s tracer]
  'vtracer, scalar, series'
  'ocean_time'
  'idVTav(itrc)'
  'v3dvar'
  1.0d0

The metadata is processed in a generic way in mod_ncparam.F. Many thanks to John Wilkin and Aboozar Tabatabai for suggesting this capability.

• Added the new C_preprocessing option LIMIT_BSTRESS to limit the bottom momentum stress in set_vbc.F to preserve bottom velocity direction, as ​suggested by John Warner. For example, for linear bottom stress we now have:

  !
  !  Set linear bottom stress.
  !
        DO j=Jstr,Jend
          DO i=IstrU,Iend
            bustr(i,j)=0.5_r8*(rdrag(i-1,j)+rdrag(i,j))*                  &
       &               u(i,j,1,nrhs)
  #    ifdef LIMIT_BSTRESS
            cff1=cff*0.5_r8*(Hz(i-1,j,1)+Hz(i,j,1))
            bustr(i,j)=SIGN(1.0_r8, bustr(i,j))*                          &
       &               MIN(ABS(bustr(i,j)),                               &
       &                   ABS(u(i,j,1,nrhs))*cff1)
  #    endif
          END DO
        END DO
        DO j=JstrV,Jend
          DO i=Istr,Iend
            bvstr(i,j)=0.5_r8*(rdrag(i,j-1)+rdrag(i,j))*                  &
       &               v(i,j,1,nrhs)
  #    ifdef LIMIT_BSTRESS
            cff1=cff*0.5_r8*(Hz(i,j-1,1)+Hz(i,j,1))
            bvstr(i,j)=SIGN(1.0_r8, bvstr(i,j))*                          &
       &               MIN(ABS(bvstr(i,j)),                               &
       &                   ABS(v(i,j,1,nrhs))*cff1)
  #    endif
          END DO
        END DO
  

  Similar strategy is used for quadratic and logarithmic formulations. I still need to work out the TLM, RPM, and ADM transformation for this routine,

• Modified get_varcoords.F to expand the processing of the coordinates attribute in input NetCDF forcing files:

          float my_var(time, lat, lon) ;
                  my_var:long_name = "my variable long name" ;
                  my_var:units = "my variable units" ;
                  my_var:coordinates = "lon lat time" ;
                  my_var:time = "time" ;
  !
  !  The following "coordinates" attribute is also allowed:
  !
                  my_var:coordinates = "lon lat" ;
  

  This is done for CF compliance.

• Corrected a bug in t3dmix4_iso.h when multiplying by the diffusion coefficient. This is a bad one! Changes are not required in the TLM, RPM, and ADM versions of this routine since these codes operate the diffusion coefficient factor at the end of the flux computation. This is done because in the linearized codes the operation needs to be recursive. Therefore, the bug in the NLM due to the split scheme is not present in the TLM, RPM, and ADM codes.

• Update all the standard input scripts in the test repository.