pkg/exf/exf_wind.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.8 2007/04/16 23:27:21 jmc Exp $
C $Name:  $

#include "EXF_OPTIONS.h"

      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"
#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)
#ifdef ALLOW_ATM_WIND
      _RL     wsSq
#else
      _RL     usSq, recip_sqrtRhoA, ustar
      _RL     tmp1, tmp2, tmp3, tmp4
      _RL     oneThirdRL
      PARAMETER ( oneThirdRL = 1.d0 / 3.d0 )
#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)

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_ATM_WIND

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

#else  /* ifndef ALLOW_ATM_WIND */

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 ; devive 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
        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 /* ifndef ALLOW_ATM_WIND */

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

C--   end bi,bj loops
       ENDDO
      ENDDO

      RETURN
      END
1	jmc	1.9	C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.8 2007/04/16 23:27:21 jmc Exp $
2	jmc	1.4	C $Name: $
3	heimbach	1.1
4			#include "EXF_OPTIONS.h"
5
6	jmc	1.9	SUBROUTINE EXF_WIND( myTime, myIter, myThid )
7	heimbach	1.1
8	jmc	1.9	C ==================================================================
9			C SUBROUTINE exf_wind
10			C ==================================================================
11			C
12			C o Prepare wind speed and stress fields
13			C
14			C ==================================================================
15			C SUBROUTINE exf_wind
16			C ==================================================================
17	heimbach	1.1
18	jmc	1.9	IMPLICIT NONE
19	heimbach	1.1
20	jmc	1.9	C == global variables ==
21	heimbach	1.1
22	jmc	1.4	#include "SIZE.h"
23	heimbach	1.1	#include "EEPARAMS.h"
24			#include "PARAMS.h"
25
26	jmc	1.8	#include "EXF_PARAM.h"
27			#include "EXF_FIELDS.h"
28			#include "EXF_CONSTANTS.h"
29	heimbach	1.1
30			#ifdef ALLOW_AUTODIFF_TAMC
31			#include "tamc.h"
32			#endif
33
34	jmc	1.9	C == routine arguments ==
35	heimbach	1.1
36	jmc	1.9	_RL myTime
37			INTEGER myIter
38			INTEGER myThid
39	heimbach	1.1
40	jmc	1.9	C == external functions ==
41	heimbach	1.1
42	jmc	1.9	C == local variables ==
43	heimbach	1.1
44	jmc	1.9	INTEGER bi,bj
45			INTEGER i,j
46			_RL wsLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47			#ifdef ALLOW_ATM_WIND
48			_RL wsSq
49			#else
50			_RL usSq, recip_sqrtRhoA, ustar
51			_RL tmp1, tmp2, tmp3, tmp4
52			_RL oneThirdRL
53			PARAMETER ( oneThirdRL = 1.d0 / 3.d0 )
54			#endif
55	heimbach	1.1
56	jmc	1.9	C == end of interface ==
57	heimbach	1.1
58	jmc	1.9	C-- Use atmospheric state to compute surface fluxes.
59	heimbach	1.1
60	jmc	1.9	C Loop over tiles.
61			DO bj = myByLo(myThid),myByHi(myThid)
62			DO bi = myBxLo(myThid),myBxHi(myThid)
63
64			C-- Initialise
65			DO j = 1,sNy
66			DO i = 1,sNx
67			wsLoc(i,j) = 0. _d 0
68	heimbach	1.1	cw(i,j,bi,bj) = 0. _d 0
69			sw(i,j,bi,bj) = 0. _d 0
70			sh(i,j,bi,bj) = 0. _d 0
71	jmc	1.9	wStress(i,j,bi,bj) = 0. _d 0
72			ENDDO
73			ENDDO
74
75	heimbach	1.1	#ifdef ALLOW_ATM_WIND
76	jmc	1.9
77			C-- Wind speed and direction.
78			DO j = 1,sNy
79			DO i = 1,sNx
80			wsSq = uwind(i,j,bi,bj)*uwind(i,j,bi,bj)
81			& + vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
82			IF ( wsSq .NE. 0. _d 0 ) THEN
83			wsLoc(i,j) = SQRT(wsSq)
84			cw(i,j,bi,bj) = uwind(i,j,bi,bj)/wsLoc(i,j)
85			sw(i,j,bi,bj) = vwind(i,j,bi,bj)/wsLoc(i,j)
86			ELSE
87			wsLoc(i,j) = 0. _d 0
88	heimbach	1.1	cw(i,j,bi,bj) = 0. _d 0
89			sw(i,j,bi,bj) = 0. _d 0
90	jmc	1.9	ENDIF
91			ENDDO
92			ENDDO
93			IF ( wspeedfile .EQ. ' ' ) THEN
94			C- wind-speed is not loaded from file: save local array into common block
95			DO j = 1,sNy
96			DO i = 1,sNx
97			wspeed(i,j,bi,bj) = wsLoc(i,j)
98			ENDDO
99			ENDDO
100			ENDIF
101
102	heimbach	1.1	#else /* ifndef ALLOW_ATM_WIND */
103	jmc	1.9
104			C-- Wind stress and direction.
105			DO j = 1,sNy
106			DO i = 1,sNx
107			IF ( stressIsOnCgrid ) THEN
108			usSq = ( ustress(i, j,bi,bj)*ustress(i ,j,bi,bj)
109			& +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
110			& +vstress(i,j, bi,bj)*vstress(i,j ,bi,bj)
111			& +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj)
112			& )*0.5 _d 0
113			ELSE
114			usSq = ustress(i,j,bi,bj)*ustress(i,j,bi,bj)
115			& +vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
116			ENDIF
117			IF ( usSq .NE. 0. _d 0 ) THEN
118			wStress(i,j,bi,bj) = SQRT(usSq)
119			c ustar = SQRT(usSq/atmrho)
120			cw(i,j,bi,bj) = ustress(i,j,bi,bj)/wStress(i,j,bi,bj)
121			sw(i,j,bi,bj) = vstress(i,j,bi,bj)/wStress(i,j,bi,bj)
122			ELSE
123			wStress(i,j,bi,bj) = 0. _d 0
124			cw(i,j,bi,bj) = 0. _d 0
125			sw(i,j,bi,bj) = 0. _d 0
126			ENDIF
127			ENDDO
128			ENDDO
129
130			IF ( wspeedfile .EQ. ' ' ) THEN
131			C-- wspeed is not loaded ; devive wind-speed by inversion of
132			C wind-stress=fct(wind-speed) relation:
133			C The variables us, sh and rdn have to be computed from
134			C given wind stresses inverting relationship for neutral
135			C drag coeff. cdn.
136			C The inversion is based on linear and quadratic form of
137			C cdn(umps); ustar can be directly computed from stress;
138			recip_sqrtRhoA = 1. _d 0 / SQRT(atmrho)
139			DO j = 1,sNy
140			DO i = 1,sNx
141			ustar = wStress(i,j,bi,bj)*recip_sqrtRhoA
142			IF ( ustar .EQ. 0. _d 0 ) THEN
143			wsLoc(i,j) = 0. _d 0
144			ELSE IF ( ustar .LT. ustofu11 ) THEN
145			tmp1 = -cquadrag_2/cquadrag_1*exf_half
146			tmp2 = SQRT(tmp1tmp1 + ustarustar/cquadrag_1)
147			wsLoc(i,j) = SQRT(tmp1 + tmp2)
148			ELSE
149			tmp1 = clindrag_2/clindrag_1*oneThirdRL
150			tmp2 = ustarustar/clindrag_1exf_half
151			& - tmp1tmp1tmp1
152			tmp3 = SQRT( ustarustar/clindrag_1
153			& (ustarustar/clindrag_10.25 _d 0 - tmp1tmp1tmp1 )
154			& )
155			tmp4 = (tmp2 + tmp3)**oneThirdRL
156			wsLoc(i,j) = tmp4 + tmp1*tmp1 / tmp4 - tmp1
157			c wsLoc(i,j) = (tmp2 + tmp3)**oneThirdRL +
158			c & tmp1tmp1 (tmp2 + tmp3)**(-oneThirdRL) - tmp1
159			ENDIF
160			ENDDO
161			ENDDO
162			C- save local array wind-speed to common block
163			DO j = 1,sNy
164			DO i = 1,sNx
165			wspeed(i,j,bi,bj) = wsLoc(i,j)
166			ENDDO
167			ENDDO
168			ENDIF
169
170			C-- infer wind field from wind-speed & wind-stress direction
171			DO j = 1,sNy
172			DO i = 1,sNx
173			uwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*cw(i,j,bi,bj)
174			vwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*sw(i,j,bi,bj)
175			ENDDO
176			ENDDO
177	heimbach	1.1	#endif /* ifndef ALLOW_ATM_WIND */
178
179	jmc	1.9	C-- set wind-speed lower limit
180			DO j = 1,sNy
181			DO i = 1,sNx
182			sh(i,j,bi,bj) = MAX(wspeed(i,j,bi,bj),uMin)
183			ENDDO
184			ENDDO
185
186			C-- end bi,bj loops
187			ENDDO
188			ENDDO
189	heimbach	1.1
190	jmc	1.9	RETURN
191			END