pkg/exf/exf_wind.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.17 2017/01/27 17:22:55 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
1	jmc	1.18	C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.17 2017/01/27 17:22:55 jmc Exp $
2	jmc	1.4	C $Name: $
3	heimbach	1.1
4			#include "EXF_OPTIONS.h"
5	gforget	1.15	#ifdef ALLOW_CTRL
6			# include "CTRL_OPTIONS.h"
7			#endif
8			#ifdef ALLOW_AUTODIFF
9			# include "AUTODIFF_OPTIONS.h"
10			#endif
11	heimbach	1.1
12	jmc	1.9	SUBROUTINE EXF_WIND( myTime, myIter, myThid )
13	heimbach	1.1
14	jmc	1.9	C ==================================================================
15			C SUBROUTINE exf_wind
16			C ==================================================================
17			C
18			C o Prepare wind speed and stress fields
19			C
20			C ==================================================================
21			C SUBROUTINE exf_wind
22			C ==================================================================
23	heimbach	1.1
24	jmc	1.9	IMPLICIT NONE
25	heimbach	1.1
26	jmc	1.9	C == global variables ==
27	heimbach	1.1
28	jmc	1.4	#include "SIZE.h"
29	heimbach	1.1	#include "EEPARAMS.h"
30			#include "PARAMS.h"
31
32	jmc	1.8	#include "EXF_PARAM.h"
33			#include "EXF_FIELDS.h"
34			#include "EXF_CONSTANTS.h"
35	heimbach	1.1
36			#ifdef ALLOW_AUTODIFF_TAMC
37			#include "tamc.h"
38	gforget	1.11	#include "tamc_keys.h"
39	heimbach	1.1	#endif
40	jmc	1.17	#ifdef ALLOW_GENTIM2D_CONTROL
41	gforget	1.13	# include "ctrl.h"
42			# include "optim.h"
43			# include "CTRL_SIZE.h"
44			# include "CTRL_GENARR.h"
45			#endif
46	heimbach	1.1
47	jmc	1.9	C == routine arguments ==
48	heimbach	1.1
49	jmc	1.9	_RL myTime
50			INTEGER myIter
51			INTEGER myThid
52	heimbach	1.1
53	jmc	1.9	C == external functions ==
54	heimbach	1.1
55	jmc	1.9	C == local variables ==
56	heimbach	1.1
57	jmc	1.9	INTEGER bi,bj
58			INTEGER i,j
59			_RL wsLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60			_RL wsSq
61	jmc	1.17	#ifdef ALLOW_BULKFORMULAE
62	jmc	1.9	_RL usSq, recip_sqrtRhoA, ustar
63			_RL tmp1, tmp2, tmp3, tmp4
64	jmc	1.17	#endif /* ALLOW_BULKFORMULAE */
65	jmc	1.9	_RL oneThirdRL
66			PARAMETER ( oneThirdRL = 1.d0 / 3.d0 )
67	jmc	1.17	#if !(defined ALLOW_BULKFORMULAE) \|\| !(defined ALLOW_ATM_TEMP)
68			_RL wsm, tmpbulk
69			#endif
70			#ifdef ALLOW_GENTIM2D_CONTROL
71	gforget	1.13	INTEGER iarr
72			#endif
73	jmc	1.9	C == end of interface ==
74	heimbach	1.1
75	jmc	1.9	C-- Use atmospheric state to compute surface fluxes.
76	heimbach	1.1
77	jmc	1.9	C Loop over tiles.
78			DO bj = myByLo(myThid),myByHi(myThid)
79			DO bi = myBxLo(myThid),myBxHi(myThid)
80
81	gforget	1.11	#ifdef ALLOW_AUTODIFF_TAMC
82			act1 = bi - myBxLo(myThid)
83			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
84			act2 = bj - myByLo(myThid)
85			max2 = myByHi(myThid) - myByLo(myThid) + 1
86			act3 = myThid - 1
87			max3 = nTx*nTy
88			act4 = ikey_dynamics - 1
89			ikey = (act1 + 1) + act2*max1
90			& + act3max1max2
91			& + act4max1max2*max3
92			#endif /* ALLOW_AUTODIFF_TAMC */
93
94	jmc	1.9	C-- Initialise
95			DO j = 1,sNy
96			DO i = 1,sNx
97			wsLoc(i,j) = 0. _d 0
98	heimbach	1.1	cw(i,j,bi,bj) = 0. _d 0
99			sw(i,j,bi,bj) = 0. _d 0
100			sh(i,j,bi,bj) = 0. _d 0
101	jmc	1.9	wStress(i,j,bi,bj) = 0. _d 0
102			ENDDO
103			ENDDO
104
105	gforget	1.12	#ifdef ALLOW_AUTODIFF_TAMC
106			CADJ STORE uwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
107			CADJ STORE vwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
108			CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
109			#endif
110
111			IF ( useAtmWind ) THEN
112	jmc	1.9
113			C-- Wind speed and direction.
114	jmc	1.17	DO j = 1,sNy
115			DO i = 1,sNx
116	jmc	1.9	wsSq = uwind(i,j,bi,bj)*uwind(i,j,bi,bj)
117	jmc	1.17	& + vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
118	jmc	1.9	IF ( wsSq .NE. 0. _d 0 ) THEN
119			wsLoc(i,j) = SQRT(wsSq)
120			cw(i,j,bi,bj) = uwind(i,j,bi,bj)/wsLoc(i,j)
121			sw(i,j,bi,bj) = vwind(i,j,bi,bj)/wsLoc(i,j)
122			ELSE
123			wsLoc(i,j) = 0. _d 0
124	heimbach	1.1	cw(i,j,bi,bj) = 0. _d 0
125			sw(i,j,bi,bj) = 0. _d 0
126	jmc	1.9	ENDIF
127	jmc	1.17	ENDDO
128	jmc	1.9	ENDDO
129	jmc	1.17	IF ( wspeedfile .EQ. ' ' ) THEN
130	jmc	1.9	C- wind-speed is not loaded from file: save local array into common block
131			DO j = 1,sNy
132			DO i = 1,sNx
133			wspeed(i,j,bi,bj) = wsLoc(i,j)
134			ENDDO
135			ENDDO
136	jmc	1.17	ENDIF
137	jmc	1.9
138	jmc	1.17	#ifdef ALLOW_BULKFORMULAE
139	gforget	1.12	ELSE
140	jmc	1.17	C-- case useAtmWind=F
141	jmc	1.9
142			C-- Wind stress and direction.
143	jmc	1.17	DO j = 1,sNy
144			DO i = 1,sNx
145	jmc	1.9	IF ( stressIsOnCgrid ) THEN
146			usSq = ( ustress(i, j,bi,bj)*ustress(i ,j,bi,bj)
147			& +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
148			& +vstress(i,j, bi,bj)*vstress(i,j ,bi,bj)
149			& +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj)
150			& )*0.5 _d 0
151			ELSE
152			usSq = ustress(i,j,bi,bj)*ustress(i,j,bi,bj)
153			& +vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
154			ENDIF
155			IF ( usSq .NE. 0. _d 0 ) THEN
156			wStress(i,j,bi,bj) = SQRT(usSq)
157			c ustar = SQRT(usSq/atmrho)
158			cw(i,j,bi,bj) = ustress(i,j,bi,bj)/wStress(i,j,bi,bj)
159			sw(i,j,bi,bj) = vstress(i,j,bi,bj)/wStress(i,j,bi,bj)
160			ELSE
161			wStress(i,j,bi,bj) = 0. _d 0
162			cw(i,j,bi,bj) = 0. _d 0
163			sw(i,j,bi,bj) = 0. _d 0
164			ENDIF
165	jmc	1.17	ENDDO
166	jmc	1.9	ENDDO
167
168	jmc	1.17	IF ( wspeedfile .EQ. ' ' ) THEN
169	mlosch	1.10	C-- wspeed is not loaded ; derive wind-speed by inversion of
170	jmc	1.9	C wind-stress=fct(wind-speed) relation:
171			C The variables us, sh and rdn have to be computed from
172			C given wind stresses inverting relationship for neutral
173			C drag coeff. cdn.
174			C The inversion is based on linear and quadratic form of
175			C cdn(umps); ustar can be directly computed from stress;
176	jmc	1.17	recip_sqrtRhoA = 1. _d 0 / SQRT(atmrho)
177			DO j = 1,sNy
178			DO i = 1,sNx
179	jmc	1.9	ustar = wStress(i,j,bi,bj)*recip_sqrtRhoA
180			IF ( ustar .EQ. 0. _d 0 ) THEN
181			wsLoc(i,j) = 0. _d 0
182			ELSE IF ( ustar .LT. ustofu11 ) THEN
183			tmp1 = -cquadrag_2/cquadrag_1*exf_half
184			tmp2 = SQRT(tmp1tmp1 + ustarustar/cquadrag_1)
185			wsLoc(i,j) = SQRT(tmp1 + tmp2)
186			ELSE
187			tmp1 = clindrag_2/clindrag_1*oneThirdRL
188			tmp2 = ustarustar/clindrag_1exf_half
189			& - tmp1tmp1tmp1
190			tmp3 = SQRT( ustarustar/clindrag_1
191			& (ustarustar/clindrag_10.25 _d 0 - tmp1tmp1tmp1 )
192			& )
193			tmp4 = (tmp2 + tmp3)**oneThirdRL
194			wsLoc(i,j) = tmp4 + tmp1*tmp1 / tmp4 - tmp1
195			c wsLoc(i,j) = (tmp2 + tmp3)**oneThirdRL +
196			c & tmp1tmp1 (tmp2 + tmp3)**(-oneThirdRL) - tmp1
197			ENDIF
198	jmc	1.17	ENDDO
199	jmc	1.9	ENDDO
200			C- save local array wind-speed to common block
201	jmc	1.17	DO j = 1,sNy
202			DO i = 1,sNx
203	jmc	1.9	wspeed(i,j,bi,bj) = wsLoc(i,j)
204	jmc	1.17	ENDDO
205	jmc	1.9	ENDDO
206	jmc	1.17	C- end if wspeedfile = empty
207			ENDIF
208	jmc	1.9
209	gforget	1.12	#ifdef ALLOW_AUTODIFF_TAMC
210			CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
211			CADJ STORE uwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
212			CADJ STORE vwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
213			CADJ STORE cw (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
214			CADJ STORE sw (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
215			#endif
216
217	jmc	1.9	C-- infer wind field from wind-speed & wind-stress direction
218	jmc	1.17	DO j = 1,sNy
219			DO i = 1,sNx
220	jmc	1.9	uwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*cw(i,j,bi,bj)
221			vwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*sw(i,j,bi,bj)
222	jmc	1.17	ENDDO
223	jmc	1.9	ENDDO
224	jmc	1.17
225			#endif /* ALLOW_BULKFORMULAE */
226			C-- end if/else useAtmWind
227	gforget	1.12	ENDIF
228	heimbach	1.1
229	jmc	1.17	#ifdef ALLOW_GENTIM2D_CONTROL
230	gforget	1.13	DO j = 1,sNy
231			DO i = 1,sNx
232			do iarr = 1, maxCtrlTim2D
233	gforget	1.16	if (xx_gentim2d_file(iarr)(1:9).EQ.'xx_wspeed')
234	gforget	1.13	& wspeed(i,j,bi,bj)=wspeed(i,j,bi,bj)+
235			& xx_gentim2d(i,j,bi,bj,iarr)
236			enddo
237			ENDDO
238			ENDDO
239			#endif
240
241	gforget	1.11	#ifdef ALLOW_AUTODIFF_TAMC
242			CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
243			#endif
244
245	jmc	1.9	C-- set wind-speed lower limit
246			DO j = 1,sNy
247			DO i = 1,sNx
248			sh(i,j,bi,bj) = MAX(wspeed(i,j,bi,bj),uMin)
249			ENDDO
250			ENDDO
251
252	jmc	1.17	#if !(defined ALLOW_BULKFORMULAE) \|\| !(defined ALLOW_ATM_TEMP)
253			C Note: In case ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP are defined,
254			C wind-stress is computed (if needed) within S/R EXF_BULKFORMULAE
255			IF ( useAtmWind ) THEN
256	jmc	1.18	c#ifdef ALLOW_ATM_WIND
257	jmc	1.17	C-- Computes wind-stress:
258			DO j = 1,sNy
259			DO i = 1,sNx
260			wsm = sh(i,j,bi,bj)
261			tmpbulk = exf_scal_BulkCdn
262			& * ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )
263			ustress(i,j,bi,bj) = atmrhotmpbulkwspeed(i,j,bi,bj)
264			& * uwind(i,j,bi,bj)
265			vstress(i,j,bi,bj) = atmrhotmpbulkwspeed(i,j,bi,bj)
266			& * vwind(i,j,bi,bj)
267			ENDDO
268			ENDDO
269	jmc	1.18	c#else /* ALLOW_ATM_WIND */
270			c STOP 'ABNORMAL END: S/R EXF_WIND: missing code for useAtmWind'
271			c#endif /* ALLOW_ATM_WIND */
272	jmc	1.17	C-- end if useAtmWind
273			ENDIF
274			#endif /* ndef ALLOW_BULKFORMULAE or ndef ALLOW_ATM_TEMP */
275
276	jmc	1.9	C-- end bi,bj loops
277			ENDDO
278			ENDDO
279	heimbach	1.1
280	jmc	1.9	RETURN
281			END