pkg/exf/exf_wind.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.14 2014/09/12 14:57:44 gforget 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
#if (defined (ALLOW_CTRL) && \
     defined (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
      _RL     usSq, recip_sqrtRhoA, ustar
      _RL     tmp1, tmp2, tmp3, tmp4
      _RL     oneThirdRL
      PARAMETER ( oneThirdRL = 1.d0 / 3.d0 )
#if (defined (ALLOW_CTRL) && \
     defined (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

        ELSE

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
        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

#if (defined (ALLOW_CTRL) && \
     defined (ALLOW_GENTIM2D_CONTROL))
        DO j = 1,sNy
         DO i = 1,sNx
           do iarr = 1, maxCtrlTim2D
           if (xx_gentim2d_file(iarr).EQ.'xxg_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

C--   end bi,bj loops
       ENDDO
      ENDDO

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