VERY IMPORTANT: surface heat/freshwater fluxes and coupling

[arango]

I modified the logic in bulk_flux.F and set_vbc.F when computing/loading the surface net heat flux and surface freshwater flux. During direct atmosphere coupling, we need to separate in different variables the fluxes from the atmosphere from ROMS surface tracer flux (vertical boundary condition) state variables stflux(:,:,itemp) and stflux(:,:,isalt).

From now on, the kinematic surface net heat flux is stored in new variable snhflx, and the surface freshwater flux is stored in EminusP. The reason for the change is to have robust simulations that are not affected by the coupling interval between atmosphere and ocean. Since ROMS manipulates stflux in set_vbc.F further according to various CPP options, it may give the wrong results if ROMS timestep is smaller than the coupling interval. In such cases, we need to persist the values of snhflx and EminusP for the coupling interval. Recall that we need to multiply surface freshwater flux (E-P) with ROMS surface salinity at every timestep. If the values are persisted, stflx(:,:,isalt) will increase by the surface salinity factor each time. Therefore, we need to re-initialize stflx(:,:,isalt) with EminusP at every timestep.

---

In set_vbc.F, we now have:

#  if defined BULK_FLUXES || defined ATM_COUPLING
!
!-----------------------------------------------------------------------
!  If air-sea parameterization and or atmospheric coupling, assign
!  surface net heat flux (degC m/s).  During model coupling, we need
!  to make sure that this forcing is unaltered during the coupling
!  interval when ROMS timestep size is smaller. Notice that further
!  manipulations to state variable stflx(:,:,itemp)are allowed below.
!-----------------------------------------------------------------------
!
      DO j=JstrR,JendR
        DO i=IstrR,IendR
          stflx(i,j,itemp)=snhflx(i,j)
        END DO
      END DO
#  endif

#  ifdef QCORRECTION
!
!-----------------------------------------------------------------------
!  Add in flux correction to surface net heat flux (degC m/s).
!-----------------------------------------------------------------------
!
! Add in net heat flux correction.
!
      DO j=JstrR,JendR
        DO i=IstrR,IendR
          stflx(i,j,itemp)=stflx(i,j,itemp)+                            &
     &                     dqdt(i,j)*(t(i,j,N(ng),nrhs,itemp)-sst(i,j))
        END DO
      END DO
#  endif

#  ifdef LIMIT_STFLX_COOLING
!
!-----------------------------------------------------------------------
!  If net heat flux is cooling and SST is at freezing point or below
!  then suppress further cooling. Note: stflx sign convention is that
!  positive means heating the ocean (J Wilkin).
!-----------------------------------------------------------------------
!
!  Below the surface heat flux stflx(:,:,itemp) is ZERO if cooling AND
!  the SST is cooler than the threshold.  The value is retained if
!  warming.
!
!    cff3 = 0      if SST warmer than threshold (cff1) - change nothing
!    cff3 = 1      if SST colder than threshold (cff1)
!
!    0.5*(cff2-ABS(cff2)) = 0                        if flux is warming
!                         = stflx(:,:,itemp)         if flux is cooling
!
      cff1=-2.0_r8              ! nominal SST threshold to cease cooling
      DO j=JstrR,JendR
        DO i=IstrR,IendR
          cff2=stflx(i,j,itemp)
          cff3=0.5_r8*(1.0_r8+SIGN(1.0_r8,cff1-t(i,j,N(ng),nrhs,itemp)))
          stflx(i,j,itemp)=cff2-cff3*0.5_r8*(cff2-ABS(cff2))
        END DO
      END DO
#  endif

#  ifdef SALINITY
!
!-----------------------------------------------------------------------
!  Multiply fresh water flux with surface salinity. If appropriate,
!  apply correction.
!-----------------------------------------------------------------------
!
      DO j=JstrR,JendR
        DO i=IstrR,IendR
#   if (defined BULK_FLUXES && defined EMINUSP) || defined ATM_COUPLING
          EmP=EminusP(i,j)
#   else
          EmP=stflx(i,j,isalt)
#   endif
#   if defined SCORRECTION
          stflx(i,j,isalt)=EmP*t(i,j,N(ng),nrhs,isalt)-                 &
     &                     Tnudg(isalt,ng)*Hz(i,j,N(ng))*               &
     &                     (t(i,j,N(ng),nrhs,isalt)-sss(i,j))
#   elif defined SRELAXATION
          stflx(i,j,isalt)=-Tnudg(isalt,ng)*Hz(i,j,N(ng))*              &
     &                     (t(i,j,N(ng),nrhs,isalt)-sss(i,j))
#   else
          stflx(i,j,isalt)=EmP*t(i,j,N(ng),nrhs,isalt)
#   endif
          btflx(i,j,isalt)=btflx(i,j,isalt)*t(i,j,1,nrhs,isalt)
        END DO
      END DO
#  endif

---

Similarly, in bulk_flux.F we now have:

      Hscale=1.0_r8/(rho0*Cp)
      cff=1.0_r8/rhow
      DO j=JstrR,JendR
        DO i=IstrR,IendR
          lrflx(i,j)=LRad(i,j)*Hscale
          lhflx(i,j)=-LHeat(i,j)*Hscale
          shflx(i,j)=-SHeat(i,j)*Hscale
          snhflx(i,j)=(srflx(i,j)+lrflx(i,j)+                           &
     &                 lhflx(i,j)+shflx(i,j))
# ifdef MASKING
          snhflx(i,j)=snhflx(i,j)*rmask(i,j)
# endif
# ifdef WET_DRY
          snhflx(i,j)=snhflx(i,j)*rmask_wet(i,j)
# endif
# ifdef EMINUSP
          evap(i,j)=LHeat(i,j)/Hlv(i,j)
#  ifdef MASKING
          evap(i,j)=evap(i,j)*rmask(i,j)
#  endif
#  ifdef WET_DRY
          evap(i,j)=evap(i,j)*rmask_wet(i,j)
#  endif
          EminusP(i,j)=cff*(evap(i,j)-rain(i,j))
#  ifdef MASKING
          EminusP(i,j)=EminusP(i,j)*rmask(i,j)
#  endif
#  ifdef WET_DRY
          EminusP(i,j)=EminusP(i,j)*rmask_wet(i,j)
#  endif
# endif
        END DO
      END DO

---

Changes are made to mod_forces.Ffor the new robust strategy.