C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.18 2017/01/28 18:29:11 jmc Exp $
C $Name: $
#include "EXF_OPTIONS.h"
#ifdef ALLOW_CTRL
# include "CTRL_OPTIONS.h"
#endif
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
SUBROUTINE EXF_WIND( myTime, myIter, myThid )
C ==================================================================
C SUBROUTINE exf_wind
C ==================================================================
C
C o Prepare wind speed and stress fields
C
C ==================================================================
C SUBROUTINE exf_wind
C ==================================================================
IMPLICIT NONE
C == global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EXF_PARAM.h"
#include "EXF_FIELDS.h"
#include "EXF_CONSTANTS.h"
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif
#ifdef ALLOW_GENTIM2D_CONTROL
# include "ctrl.h"
# include "optim.h"
# include "CTRL_SIZE.h"
# include "CTRL_GENARR.h"
#endif
C == routine arguments ==
_RL myTime
INTEGER myIter
INTEGER myThid
C == external functions ==
C == local variables ==
INTEGER bi,bj
INTEGER i,j
_RL wsLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL wsSq
#ifdef ALLOW_BULKFORMULAE
_RL usSq, recip_sqrtRhoA, ustar
_RL tmp1, tmp2, tmp3, tmp4
#endif /* ALLOW_BULKFORMULAE */
_RL oneThirdRL
PARAMETER ( oneThirdRL = 1.d0 / 3.d0 )
#if !(defined ALLOW_BULKFORMULAE) || !(defined ALLOW_ATM_TEMP)
_RL wsm, tmpbulk
#endif
#ifdef ALLOW_GENTIM2D_CONTROL
INTEGER iarr
#endif
C == end of interface ==
C-- Use atmospheric state to compute surface fluxes.
C Loop over tiles.
DO bj = myByLo(myThid),myByHi(myThid)
DO bi = myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
ikey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
C-- Initialise
DO j = 1,sNy
DO i = 1,sNx
wsLoc(i,j) = 0. _d 0
cw(i,j,bi,bj) = 0. _d 0
sw(i,j,bi,bj) = 0. _d 0
sh(i,j,bi,bj) = 0. _d 0
wStress(i,j,bi,bj) = 0. _d 0
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE uwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE vwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
IF ( useAtmWind ) THEN
C-- Wind speed and direction.
DO j = 1,sNy
DO i = 1,sNx
wsSq = uwind(i,j,bi,bj)*uwind(i,j,bi,bj)
& + vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
IF ( wsSq .NE. 0. _d 0 ) THEN
wsLoc(i,j) = SQRT(wsSq)
cw(i,j,bi,bj) = uwind(i,j,bi,bj)/wsLoc(i,j)
sw(i,j,bi,bj) = vwind(i,j,bi,bj)/wsLoc(i,j)
ELSE
wsLoc(i,j) = 0. _d 0
cw(i,j,bi,bj) = 0. _d 0
sw(i,j,bi,bj) = 0. _d 0
ENDIF
ENDDO
ENDDO
IF ( wspeedfile .EQ. ' ' ) THEN
C- wind-speed is not loaded from file: save local array into common block
DO j = 1,sNy
DO i = 1,sNx
wspeed(i,j,bi,bj) = wsLoc(i,j)
ENDDO
ENDDO
ENDIF
#ifdef ALLOW_BULKFORMULAE
ELSE
C-- case useAtmWind=F
C-- Wind stress and direction.
DO j = 1,sNy
DO i = 1,sNx
IF ( stressIsOnCgrid ) THEN
usSq = ( ustress(i, j,bi,bj)*ustress(i ,j,bi,bj)
& +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
& +vstress(i,j, bi,bj)*vstress(i,j ,bi,bj)
& +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj)
& )*0.5 _d 0
ELSE
usSq = ustress(i,j,bi,bj)*ustress(i,j,bi,bj)
& +vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
ENDIF
IF ( usSq .NE. 0. _d 0 ) THEN
wStress(i,j,bi,bj) = SQRT(usSq)
c ustar = SQRT(usSq/atmrho)
cw(i,j,bi,bj) = ustress(i,j,bi,bj)/wStress(i,j,bi,bj)
sw(i,j,bi,bj) = vstress(i,j,bi,bj)/wStress(i,j,bi,bj)
ELSE
wStress(i,j,bi,bj) = 0. _d 0
cw(i,j,bi,bj) = 0. _d 0
sw(i,j,bi,bj) = 0. _d 0
ENDIF
ENDDO
ENDDO
IF ( wspeedfile .EQ. ' ' ) THEN
C-- wspeed is not loaded ; derive wind-speed by inversion of
C wind-stress=fct(wind-speed) relation:
C The variables us, sh and rdn have to be computed from
C given wind stresses inverting relationship for neutral
C drag coeff. cdn.
C The inversion is based on linear and quadratic form of
C cdn(umps); ustar can be directly computed from stress;
recip_sqrtRhoA = 1. _d 0 / SQRT(atmrho)
DO j = 1,sNy
DO i = 1,sNx
ustar = wStress(i,j,bi,bj)*recip_sqrtRhoA
IF ( ustar .EQ. 0. _d 0 ) THEN
wsLoc(i,j) = 0. _d 0
ELSE IF ( ustar .LT. ustofu11 ) THEN
tmp1 = -cquadrag_2/cquadrag_1*exf_half
tmp2 = SQRT(tmp1*tmp1 + ustar*ustar/cquadrag_1)
wsLoc(i,j) = SQRT(tmp1 + tmp2)
ELSE
tmp1 = clindrag_2/clindrag_1*oneThirdRL
tmp2 = ustar*ustar/clindrag_1*exf_half
& - tmp1*tmp1*tmp1
tmp3 = SQRT( ustar*ustar/clindrag_1*
& (ustar*ustar/clindrag_1*0.25 _d 0 - tmp1*tmp1*tmp1 )
& )
tmp4 = (tmp2 + tmp3)**oneThirdRL
wsLoc(i,j) = tmp4 + tmp1*tmp1 / tmp4 - tmp1
c wsLoc(i,j) = (tmp2 + tmp3)**oneThirdRL +
c & tmp1*tmp1 * (tmp2 + tmp3)**(-oneThirdRL) - tmp1
ENDIF
ENDDO
ENDDO
C- save local array wind-speed to common block
DO j = 1,sNy
DO i = 1,sNx
wspeed(i,j,bi,bj) = wsLoc(i,j)
ENDDO
ENDDO
C- end if wspeedfile = empty
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE uwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE vwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE cw (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE sw (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
C-- infer wind field from wind-speed & wind-stress direction
DO j = 1,sNy
DO i = 1,sNx
uwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*cw(i,j,bi,bj)
vwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*sw(i,j,bi,bj)
ENDDO
ENDDO
#endif /* ALLOW_BULKFORMULAE */
C-- end if/else useAtmWind
ENDIF
#ifdef ALLOW_GENTIM2D_CONTROL
DO j = 1,sNy
DO i = 1,sNx
do iarr = 1, maxCtrlTim2D
if (xx_gentim2d_file(iarr)(1:9).EQ.'xx_wspeed')
& wspeed(i,j,bi,bj)=wspeed(i,j,bi,bj)+
& xx_gentim2d(i,j,bi,bj,iarr)
enddo
ENDDO
ENDDO
#endif
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
C-- set wind-speed lower limit
DO j = 1,sNy
DO i = 1,sNx
sh(i,j,bi,bj) = MAX(wspeed(i,j,bi,bj),uMin)
ENDDO
ENDDO
#if !(defined ALLOW_BULKFORMULAE) || !(defined ALLOW_ATM_TEMP)
C Note: In case ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP are defined,
C wind-stress is computed (if needed) within S/R EXF_BULKFORMULAE
IF ( useAtmWind ) THEN
c#ifdef ALLOW_ATM_WIND
C-- Computes wind-stress:
DO j = 1,sNy
DO i = 1,sNx
wsm = sh(i,j,bi,bj)
tmpbulk = exf_scal_BulkCdn
& * ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )
ustress(i,j,bi,bj) = atmrho*tmpbulk*wspeed(i,j,bi,bj)
& * uwind(i,j,bi,bj)
vstress(i,j,bi,bj) = atmrho*tmpbulk*wspeed(i,j,bi,bj)
& * vwind(i,j,bi,bj)
ENDDO
ENDDO
c#else /* ALLOW_ATM_WIND */
c STOP 'ABNORMAL END: S/R EXF_WIND: missing code for useAtmWind'
c#endif /* ALLOW_ATM_WIND */
C-- end if useAtmWind
ENDIF
#endif /* ndef ALLOW_BULKFORMULAE or ndef ALLOW_ATM_TEMP */
C-- end bi,bj loops
ENDDO
ENDDO
RETURN
END