/[MITgcm]/MITgcm/pkg/exf/exf_wind.F
ViewVC logotype

Contents of /MITgcm/pkg/exf/exf_wind.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.12 - (show annotations) (download)
Tue Aug 28 19:17:46 2012 UTC (11 years, 9 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63s, checkpoint64, checkpoint65, checkpoint65b, checkpoint65c, checkpoint65a
Changes since 1.11: +18 -7 lines 
- pkg/exf : added run time switch useAtmWind to replace ALLOW_ATM_WIND
  cpp switch. ALLOW_ATM_WIND now just sets the useAtmWind default (see
  exf_readparms.F) and force defines ALLOW_BULKFORMULAE (EXF_OPTIONS.h).
- pkg/exf, autodiff, ctrl, ecco and seaice : remove ALLOW_ATM_WIND
  brackets, or replace them with useAtmWind ones.
- pkg/ctrl, ecco : allow to compile both ALLOW_U/VSTRESS_CONTROL and
  ALLOW_U/VWIND_CONTROL. Depending on useAtmWind, one is inactive,
  and the other is active (see exf_getffields.F/exf_getsurfacefluxes.F).
1 C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_wind.F,v 1.11 2010/10/25 16:21:58 gforget Exp $
2 C $Name: $
3
4 #include "EXF_OPTIONS.h"
5
6 SUBROUTINE EXF_WIND( myTime, myIter, myThid )
7
8 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
18 IMPLICIT NONE
19
20 C == global variables ==
21
22 #include "SIZE.h"
23 #include "EEPARAMS.h"
24 #include "PARAMS.h"
25
26 #include "EXF_PARAM.h"
27 #include "EXF_FIELDS.h"
28 #include "EXF_CONSTANTS.h"
29
30 #ifdef ALLOW_AUTODIFF_TAMC
31 #include "tamc.h"
32 #include "tamc_keys.h"
33 #endif
34
35 C == routine arguments ==
36
37 _RL myTime
38 INTEGER myIter
39 INTEGER myThid
40
41 C == external functions ==
42
43 C == local variables ==
44
45 INTEGER bi,bj
46 INTEGER i,j
47 _RL wsLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48 _RL wsSq
49 _RL usSq, recip_sqrtRhoA, ustar
50 _RL tmp1, tmp2, tmp3, tmp4
51 _RL oneThirdRL
52 PARAMETER ( oneThirdRL = 1.d0 / 3.d0 )
53
54 C == end of interface ==
55
56 C-- Use atmospheric state to compute surface fluxes.
57
58 C Loop over tiles.
59 DO bj = myByLo(myThid),myByHi(myThid)
60 DO bi = myBxLo(myThid),myBxHi(myThid)
61
62 #ifdef ALLOW_AUTODIFF_TAMC
63 act1 = bi - myBxLo(myThid)
64 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
65 act2 = bj - myByLo(myThid)
66 max2 = myByHi(myThid) - myByLo(myThid) + 1
67 act3 = myThid - 1
68 max3 = nTx*nTy
69 act4 = ikey_dynamics - 1
70 ikey = (act1 + 1) + act2*max1
71 & + act3*max1*max2
72 & + act4*max1*max2*max3
73 #endif /* ALLOW_AUTODIFF_TAMC */
74
75 C-- Initialise
76 DO j = 1,sNy
77 DO i = 1,sNx
78 wsLoc(i,j) = 0. _d 0
79 cw(i,j,bi,bj) = 0. _d 0
80 sw(i,j,bi,bj) = 0. _d 0
81 sh(i,j,bi,bj) = 0. _d 0
82 wStress(i,j,bi,bj) = 0. _d 0
83 ENDDO
84 ENDDO
85
86 #ifdef ALLOW_AUTODIFF_TAMC
87 CADJ STORE uwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
88 CADJ STORE vwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
89 CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
90 #endif
91
92 IF ( useAtmWind ) THEN
93
94 C-- Wind speed and direction.
95 DO j = 1,sNy
96 DO i = 1,sNx
97 wsSq = uwind(i,j,bi,bj)*uwind(i,j,bi,bj)
98 & + vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
99 IF ( wsSq .NE. 0. _d 0 ) THEN
100 wsLoc(i,j) = SQRT(wsSq)
101 cw(i,j,bi,bj) = uwind(i,j,bi,bj)/wsLoc(i,j)
102 sw(i,j,bi,bj) = vwind(i,j,bi,bj)/wsLoc(i,j)
103 ELSE
104 wsLoc(i,j) = 0. _d 0
105 cw(i,j,bi,bj) = 0. _d 0
106 sw(i,j,bi,bj) = 0. _d 0
107 ENDIF
108 ENDDO
109 ENDDO
110 IF ( wspeedfile .EQ. ' ' ) THEN
111 C- wind-speed is not loaded from file: save local array into common block
112 DO j = 1,sNy
113 DO i = 1,sNx
114 wspeed(i,j,bi,bj) = wsLoc(i,j)
115 ENDDO
116 ENDDO
117 ENDIF
118
119 ELSE
120
121 C-- Wind stress and direction.
122 DO j = 1,sNy
123 DO i = 1,sNx
124 IF ( stressIsOnCgrid ) THEN
125 usSq = ( ustress(i, j,bi,bj)*ustress(i ,j,bi,bj)
126 & +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
127 & +vstress(i,j, bi,bj)*vstress(i,j ,bi,bj)
128 & +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj)
129 & )*0.5 _d 0
130 ELSE
131 usSq = ustress(i,j,bi,bj)*ustress(i,j,bi,bj)
132 & +vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
133 ENDIF
134 IF ( usSq .NE. 0. _d 0 ) THEN
135 wStress(i,j,bi,bj) = SQRT(usSq)
136 c ustar = SQRT(usSq/atmrho)
137 cw(i,j,bi,bj) = ustress(i,j,bi,bj)/wStress(i,j,bi,bj)
138 sw(i,j,bi,bj) = vstress(i,j,bi,bj)/wStress(i,j,bi,bj)
139 ELSE
140 wStress(i,j,bi,bj) = 0. _d 0
141 cw(i,j,bi,bj) = 0. _d 0
142 sw(i,j,bi,bj) = 0. _d 0
143 ENDIF
144 ENDDO
145 ENDDO
146
147 IF ( wspeedfile .EQ. ' ' ) THEN
148 C-- wspeed is not loaded ; derive wind-speed by inversion of
149 C wind-stress=fct(wind-speed) relation:
150 C The variables us, sh and rdn have to be computed from
151 C given wind stresses inverting relationship for neutral
152 C drag coeff. cdn.
153 C The inversion is based on linear and quadratic form of
154 C cdn(umps); ustar can be directly computed from stress;
155 recip_sqrtRhoA = 1. _d 0 / SQRT(atmrho)
156 DO j = 1,sNy
157 DO i = 1,sNx
158 ustar = wStress(i,j,bi,bj)*recip_sqrtRhoA
159 IF ( ustar .EQ. 0. _d 0 ) THEN
160 wsLoc(i,j) = 0. _d 0
161 ELSE IF ( ustar .LT. ustofu11 ) THEN
162 tmp1 = -cquadrag_2/cquadrag_1*exf_half
163 tmp2 = SQRT(tmp1*tmp1 + ustar*ustar/cquadrag_1)
164 wsLoc(i,j) = SQRT(tmp1 + tmp2)
165 ELSE
166 tmp1 = clindrag_2/clindrag_1*oneThirdRL
167 tmp2 = ustar*ustar/clindrag_1*exf_half
168 & - tmp1*tmp1*tmp1
169 tmp3 = SQRT( ustar*ustar/clindrag_1*
170 & (ustar*ustar/clindrag_1*0.25 _d 0 - tmp1*tmp1*tmp1 )
171 & )
172 tmp4 = (tmp2 + tmp3)**oneThirdRL
173 wsLoc(i,j) = tmp4 + tmp1*tmp1 / tmp4 - tmp1
174 c wsLoc(i,j) = (tmp2 + tmp3)**oneThirdRL +
175 c & tmp1*tmp1 * (tmp2 + tmp3)**(-oneThirdRL) - tmp1
176 ENDIF
177 ENDDO
178 ENDDO
179 C- save local array wind-speed to common block
180 DO j = 1,sNy
181 DO i = 1,sNx
182 wspeed(i,j,bi,bj) = wsLoc(i,j)
183 ENDDO
184 ENDDO
185 ENDIF
186
187 #ifdef ALLOW_AUTODIFF_TAMC
188 CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
189 CADJ STORE uwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
190 CADJ STORE vwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
191 CADJ STORE cw (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
192 CADJ STORE sw (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
193 #endif
194
195 C-- infer wind field from wind-speed & wind-stress direction
196 DO j = 1,sNy
197 DO i = 1,sNx
198 uwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*cw(i,j,bi,bj)
199 vwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*sw(i,j,bi,bj)
200 ENDDO
201 ENDDO
202 ENDIF
203
204 #ifdef ALLOW_AUTODIFF_TAMC
205 CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
206 #endif
207
208 C-- set wind-speed lower limit
209 DO j = 1,sNy
210 DO i = 1,sNx
211 sh(i,j,bi,bj) = MAX(wspeed(i,j,bi,bj),uMin)
212 ENDDO
213 ENDDO
214
215 C-- end bi,bj loops
216 ENDDO
217 ENDDO
218
219 RETURN
220 END
  ViewVC Help
Powered by ViewVC 1.1.22