/[MITgcm]/MITgcm_contrib/verification_other/offline_cheapaml/code/cheapaml.F
ViewVC logotype

Contents of /MITgcm_contrib/verification_other/offline_cheapaml/code/cheapaml.F

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

Revision 1.2 - (show annotations) (download)
Sat May 25 18:19:05 2013 UTC (12 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64h
Changes since 1.1: +119 -64 lines 
- consistent treatment of uWind,vWind location: assume everywhere in
  pkg/cheapaml that they are on C-grid (@ uVel,vVel location)
  (Note: was already the case for Tair,Qair advection);
- fill in array "Qsw" (for short-wave heating);
- add diagnostic for precip ;
1 C $Header: /u/gcmpack/MITgcm_contrib/verification_other/offline_cheapaml/code/cheapaml.F,v 1.1 2013/05/22 19:39:51 jmc Exp $
2 C $Name: $
3
4 #include "CHEAPAML_OPTIONS.h"
5
6 SUBROUTINE CHEAPAML(
7 I myTime, myIter, myThid )
8
9 C ==================================================================
10 C SUBROUTINE cheapaml
11 C ==================================================================
12 C
13 C o Get the surface fluxes used to force ocean model
14 C
15 C Output:
16 C ------
17 C ustress, vstress - wind stress
18 C Qnet - net heat flux
19 C EmPmR - net freshwater flux
20 C Tair - mean air temperature (K) at height ht (m)
21 C Qair - Specific humidity kg/kg
22 C Cheaptracer - passive tracer
23 C ---------
24 C
25 C Input:
26 C ------
27 C uWind, vWind - mean wind speed (m/s)
28 C Tr - Relaxation profile for Tair on boundaries (C)
29 C qr - Relaxation profile for specific humidity (kg/kg)
30 C CheaptracerR - Relaxation profile for passive tracer
31 C ==================================================================
32 C SUBROUTINE cheapaml
33 C ==================================================================
34
35 IMPLICIT NONE
36
37 C == global variables ==
38 #include "SIZE.h"
39 #include "EEPARAMS.h"
40 #include "EESUPPORT.h"
41 #include "PARAMS.h"
42 #include "DYNVARS.h"
43 #include "GRID.h"
44 #include "FFIELDS.h"
45 #include "CHEAPAML.h"
46
47 C == routine arguments ==
48 _RL myTime
49 INTEGER myIter
50 INTEGER myThid
51
52 C == Local variables ==
53 INTEGER bi,bj
54 INTEGER i,j, nt, startAB
55 INTEGER iMin, iMax
56 INTEGER jMin, jMax
57 LOGICAL writeDbug
58 CHARACTER*10 sufx
59 LOGICAL xIsPeriodic, yIsPeriodic
60
61 C tendencies of atmospheric temperature, current and past
62 _RL gTair(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
63 _RL gqair(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
64 _RL gCheaptracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
65 C zonal and meridional transports
66 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67 _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 C AML timestep
69 _RL deltaTTracer,deltaTm,ts,xalwu
70 _RL dm,pt,xalwd,xlwnet
71 _RL dtemp,xflu,xfld,dq,dtr
72 c _RL Fclouds, ttt2
73 _RL q,precip,ttt,entrain
74 _RL uRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 _RL vRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 _RL windSq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77 _RL fsha(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78 _RL flha(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79 _RL evp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80 _RL xolw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81 _RL ssqt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82 _RL q100(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83 _RL cdq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84 C surfDrag :: surface drag coeff (for wind stress)
85 _RL surfDrag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
86 _RL dumArg(6)
87 _RL fsha0, flha0, evp_0, xolw0, ssqt0, q100_0, cdq_0
88 _RL Tsurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL iceFrac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 _RL icFrac, opFrac
91 C temp var
92 _RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93
94 C variables for htflux
95 #ifdef ALLOW_SEAGER
96 integer iperx
97 integer lsm(snx,sny)
98
99 real slat, salt_fixed,tstep
100 real dyd_htf(sny),dxd_htf(snx,sny)
101 real sst_htf(snx,sny)
102 real cldfr_htf(snx,sny),wspd_htf(snx,sny),u_htf(snx,sny)
103 real v_htf(snx,sny)
104 real q_htf(snx,sny),t_htf(snx,sny),rlh(snx,sny)
105 real sh(snx,sny),qlw(snx,sny)
106 real qsw_htf(snx,sny),ppo(snx,sny),qa(snx,sny),th(snx,sny)
107 real rh(snx,sny)
108 real qisw(snx,sny),ppi(snx,sny),hice(snx,sny),cice(snx,sny)
109 real thice(snx,sny),tsnw(snx,sny),qios(snx,sny),brne(snx,sny)
110 real rlhi(snx,sny),shi(snx,sny),qlwi(snx,sny),qswi(snx,sny)
111 real albedo(snx,sny)
112 _RL SH_sauv(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
113 _RL LH_sauv(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
114 #endif /* ALLOW_SEAGER */
115
116 C useful values
117 C inverse of time step
118 deltaTm=1. _d 0/deltaT
119
120 C atmospheric timestep
121 deltaTtracer = deltaT/FLOAT(cheapaml_ntim)
122
123 c writeDbug = debugLevel.GE.debLevC .AND.
124 c & DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock)
125 writeDbug = debugLevel.GE.debLevC .AND. diagFreq.GT.0.
126
127 #ifdef ALLOW_DIAGNOSTICS
128 C-- fill-in diagnostics for cheapAML state variables:
129 IF ( useDiagnostics ) THEN
130 CALL DIAGNOSTICS_FILL( Tair, 'CH_TAIR ', 0,1,0,1,1, myThid )
131 IF (useFreshWaterFlux)
132 & CALL DIAGNOSTICS_FILL( Qair, 'CH_QAIR ', 0,1,0,1,1, myThid )
133 ENDIF
134 #endif /* ALLOW_DIAGNOSTICS */
135
136 #ifdef ALLOW_SEAGER
137
138 C initialize array for the seager computation
139 slat = ygOrigin
140 salt_fixed = 35.0
141 iperx = 0
142 tstep = deltaT
143
144 DO bj=myByLo(myThid),myByHi(myThid)
145 DO bi=myBxLo(myThid),myBxHi(myThid)
146 DO j = 1,sny
147 DO i = 1,snx
148 C inputs
149 lsm (i,j) = 1-maskC(i,j,1,bi,bj)
150 lsm(1,j) = 1.0
151 lsm(snx,j) = 1.0
152 lsm(i,1) = 1.0
153 lsm(i,sny) = 1.0
154 c if (i.le.100) lsm(i,j) = 1.0
155
156 dyd_htf(j) = delY(j)
157 dxd_htf(i,j) = delX(i)
158 sst_htf(i,j) = theta(i,j,1,bi,bj) + celsius2K
159 cldfr_htf(i,j) = 0.0 _d 0
160 u_htf(i,j) = uwind(i,j,bi,bj)
161 v_htf(i,j) = vwind(i,j,bi,bj)
162 q_htf(i,j) = qair(i,j,bi,bj)
163 t_htf(i,j) = Tair(i,j,bi,bj) + celsius2K
164 qisw(i,j) = solar(i,j,bi,bj)
165 ppi(i,j) = 0.0 _d 0
166 wspd_htf(i,j) = sqrt(uwind(i,j,bi,bj)**2
167 & + vwind(i,j,bi,bj)**2)
168
169 cice(i,j) = 0.0 _d 0
170 C je met la temperature de la glace la dedans
171 tsnw(i,j) = 0.0 _d 0 + celsius2K
172
173 C outputs
174 C rlh(snx,sny)
175 C sh(snx,sny)
176 C qlw(snx,sny)
177 C qsw_htf(snx,sny)
178 C ppo(snx,sny)
179 C qa(snx,sny)
180 C th(snx,sny)
181 C rh(snx,sny)
182 C hice(snx,sny)
183 C thice(snx,sny)
184 C qios(snx,sny)
185 C brne(snx,sny)
186 C rlhi(snx,sny)
187 C shi(snx,sny)
188 C qlwi(snx,sny)
189 C qswi(snx,sny)
190 C albedo(snx,sny) = 0. _d 0
191 ENDDO
192 ENDDO
193
194 C close bi, bj loops
195 ENDDO
196 ENDDO
197
198 CALL HTFLUX
199 call htfluxice(snx,sny,lsm,dxd_htf,dyd_htf,tstep,
200 + sst_htf,cldfr_htf,wspd_htf,u_htf,v_htf,q_htf,t_htf
201 $ ,rlh,sh,qlw,qsw_htf,ppo,qa,th,rh,
202 + qisw,ppi,hice,cice,thice,tsnw,qios,brne,rlhi,shi,qlwi,qswi,
203 + iperx,salt_fixed,albedo,slat)
204
205 DO bj=myByLo(myThid),myByHi(myThid)
206 DO bi=myBxLo(myThid),myBxHi(myThid)
207 DO j = 1,sny
208 DO i = 1,snx
209 C OUTPUT
210 if (lsm(i,j).eq.0) then
211 qair(i,j,bi,bj) = qa(i,j)
212 Tair(i,j,bi,bj) = th(i,j) - celsius2K
213 SH_sauv(i,j,bi,bj) = sh(i,j)
214 LH_sauv(i,j,bi,bj) = rlh(i,j)
215 else
216 qair(i,j,bi,bj) = qr(i,j,bi,bj)
217 Tair(i,j,bi,bj) = tr(i,j,bi,bj)
218 endif
219
220 ENDDO
221 ENDDO
222
223 C close bi, bj loops
224 ENDDO
225 ENDDO
226
227 #else /* ALLOW_SEAGER */
228
229 DO bj=myByLo(myThid),myByHi(myThid)
230 DO bi=myBxLo(myThid),myBxHi(myThid)
231 C initialize net heat flux and fresh water flux arrays
232 DO j = 1-OLy,sNy+OLy
233 DO i = 1-OLx,sNx+OLx
234 Qnet(i,j,bi,bj) = 0. _d 0
235 EmPmR(i,j,bi,bj)= 0. _d 0
236 ENDDO
237 ENDDO
238 ENDDO
239 ENDDO
240
241 C this is a reprogramming to speed up cheapaml
242 C the short atmospheric time step is applied to
243 C advection and diffusion only. diabatic forcing is computed
244 C once and used for the entire oceanic time step.
245
246 C cycle through atmospheric advective/diffusive
247 C surface temperature evolution
248
249 DO nt=1,cheapaml_ntim
250
251 DO bj=myByLo(myThid),myByHi(myThid)
252 DO bi=myBxLo(myThid),myBxHi(myThid)
253
254 C compute advective and diffusive flux divergence
255 DO j=1-OLy,sNy+OLy
256 DO i=1-OLx,sNx+OLx
257 gTair(i,j,bi,bj)=0. _d 0
258 uTrans(i,j)=uWind(i,j,bi,bj)
259 vTrans(i,j)=vWind(i,j,bi,bj)
260 ENDDO
261 ENDDO
262 CALL GAD_2d_CALC_RHS(
263 I bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy,
264 I uTrans, vTrans,
265 I uWind, vWind,
266 I cheapaml_kdiff, Tair,
267 I deltaTtracer, zu, useFluxLimit,
268 I cheapamlXperiodic, cheapamlYperiodic,
269 O wWind,
270 U gTair,
271 I myTime, myIter, myThid )
272 startAB = cheapTairStartAB + nt - 1
273 CALL ADAMS_BASHFORTH2(
274 I bi, bj, 1, 1,
275 U gTair, gTairm, tmpFld,
276 I startAB, myIter, myThid )
277 CALL CHEAPAML_TIMESTEP(
278 I bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy, deltaTtracer,
279 I gTair,
280 U Tair,
281 I nt, myIter, myThid )
282 C close bi,bj loops
283 ENDDO
284 ENDDO
285 C update edges
286 _EXCH_XY_RL(Tair,myThid)
287
288 IF (useFreshWaterFlux) THEN
289 C do water
290 DO bj=myByLo(myThid),myByHi(myThid)
291 DO bi=myBxLo(myThid),myBxHi(myThid)
292 DO j=1-OLy,sNy+OLy
293 DO i=1-OLx,sNx+OLx
294 gqair(i,j,bi,bj)=0. _d 0
295 uTrans(i,j)=uWind(i,j,bi,bj)
296 vTrans(i,j)=vWind(i,j,bi,bj)
297 ENDDO
298 ENDDO
299 CALL GAD_2d_CALC_RHS(
300 I bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy,
301 I uTrans, vTrans,
302 I uWind, vWind,
303 I cheapaml_kdiff, qair,
304 I deltaTtracer, zu, useFluxLimit,
305 I cheapamlXperiodic, cheapamlYperiodic,
306 O wWind,
307 U gqair,
308 I myTime, myIter, myThid )
309 startAB = cheapTairStartAB + nt - 1
310 CALL ADAMS_BASHFORTH2(
311 I bi, bj, 1, 1,
312 U gqair, gqairm, tmpFld,
313 I startAB, myIter, myThid )
314 CALL CHEAPAML_TIMESTEP(
315 I bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy, deltaTtracer,
316 I gqair,
317 U qair,
318 I nt, myIter, myThid )
319 C close bi, bj loops
320 ENDDO
321 ENDDO
322 C update edges
323 _EXCH_XY_RL(qair,myThid)
324 ENDIF ! if use freshwater
325
326 IF (useCheapTracer) THEN
327 C do tracer
328 DO bj=myByLo(myThid),myByHi(myThid)
329 DO bi=myBxLo(myThid),myBxHi(myThid)
330 DO j=1-OLy,sNy+OLy
331 DO i=1-OLx,sNx+OLx
332 gCheaptracer(i,j,bi,bj)=0. _d 0
333 uTrans(i,j)=uWind(i,j,bi,bj)
334 vTrans(i,j)=vWind(i,j,bi,bj)
335 ENDDO
336 ENDDO
337 CALL GAD_2d_CALC_RHS(
338 I bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy,
339 I uTrans, vTrans,
340 I uWind, vWind,
341 I cheapaml_kdiff, Cheaptracer,
342 I deltaTtracer, zu, useFluxLimit,
343 I cheapamlXperiodic, cheapamlYperiodic,
344 O wWind,
345 U gCheaptracer,
346 I myTime, myIter, myThid )
347 startAB = cheapTracStartAB + nt - 1
348 CALL ADAMS_BASHFORTH2(
349 I bi, bj, 1, 1,
350 U gCheaptracer, gCheaptracerm, tmpFld,
351 I startAB, myIter, myThid )
352 CALL CHEAPAML_TIMESTEP(
353 I bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy, deltaTtracer,
354 I gCheaptracer,
355 U Cheaptracer,
356 I nt, myIter, myThid )
357 C close bi, bj loops
358 ENDDO
359 ENDDO
360 C update edges
361 _EXCH_XY_RL(Cheaptracer,myThid)
362 ENDIF ! if use tracer
363
364 C reset boundaries to open boundary profile
365 IF ( .NOT.(cheapamlXperiodic.AND.cheapamlYperiodic) ) THEN
366 DO bj=myByLo(myThid),myByHi(myThid)
367 DO bi=myBxLo(myThid),myBxHi(myThid)
368 CALL CHEAPAML_COPY_EDGES(
369 I cheapamlXperiodic, cheapamlYperiodic,
370 I Tr(1-OLx,1-OLy,bi,bj),
371 U Tair(1-OLx,1-OLy,bi,bj),
372 I bi, bj, myIter, myThid )
373 IF (useFreshWaterFlux) THEN
374 CALL CHEAPAML_COPY_EDGES(
375 I cheapamlXperiodic, cheapamlYperiodic,
376 I qr(1-OLx,1-OLy,bi,bj),
377 U qair(1-OLx,1-OLy,bi,bj),
378 I bi, bj, myIter, myThid )
379 ENDIF
380 IF (useCheapTracer) THEN
381 CALL CHEAPAML_COPY_EDGES(
382 I cheapamlXperiodic, cheapamlYperiodic,
383 I CheaptracerR(1-OLx,1-OLy,bi,bj),
384 U Cheaptracer(1-OLx,1-OLy,bi,bj),
385 I bi, bj, myIter, myThid )
386 ENDIF
387 ENDDO
388 ENDDO
389 ENDIF
390
391 C-- end loop on nt (short time-step loop)
392 ENDDO
393 IF ( writeDbug ) THEN
394 WRITE(sufx,'(I10.10)') myIter
395 CALL WRITE_FLD_XY_RL('tAir_afAdv.', sufx, Tair, myIter, myThid )
396 IF (useFreshWaterFlux)
397 & CALL WRITE_FLD_XY_RL('qAir_afAdv.', sufx, qair, myIter, myThid )
398 ENDIF
399
400 C cycling on short atmospheric time step is now done
401
402 C now continue with diabatic forcing
403 iMin = 1
404 iMax = sNx
405 jMin = 1
406 jMax = sNy
407
408 DO bj=myByLo(myThid),myByHi(myThid)
409 DO bi=myBxLo(myThid),myBxHi(myThid)
410
411 DO j = 1-OLy, sNy+OLy
412 DO i = 1-OLx, sNx+OLx
413 surfDrag(i,j,bi,bj) = 0.
414 uRelWind(i,j) = uWind(i,j,bi,bj)-uVel(i,j,1,bi,bj)
415 vRelWind(i,j) = vWind(i,j,bi,bj)-vVel(i,j,1,bi,bj)
416 ENDDO
417 ENDDO
418 DO j = jMin,jMax
419 DO i = iMin,iMax
420 windSq(i,j) = ( uRelWind( i ,j)*uRelWind( i ,j)
421 & + uRelWind(i+1,j)*uRelWind(i+1,j)
422 & + vRelWind(i, j )*vRelWind(i, j )
423 & + vRelWind(i,j+1)*vRelWind(i,j+1)
424 & )*0.5 _d 0
425 #ifdef INCONSISTENT_WIND_LOCATION
426 windSq(i,j) = uRelWind(i,j)*uRelWind(i,j)
427 & + vRelWind(i,j)*vRelWind(i,j)
428 #endif
429 ENDDO
430 ENDDO
431
432 IF ( useThSIce ) THEN
433 CALL CHEAPAML_SEAICE(
434 I solar(1-OLx,1-OLy,bi,bj),
435 I cheapdlongwave(1-OLx,1-OLy,bi,bj),
436 I uWind(1-OLx,1-OLy,bi,bj),
437 I vWind(1-OLx,1-OLy,bi,bj), lath,
438 O fsha, flha, evp, xolw, ssqt, q100, cdq,
439 O Tsurf, iceFrac, Qsw(1-OLx,1-OLy,bi,bj),
440 I bi, bj, myTime, myIter, myThid )
441 DO j = jMin,jMax
442 DO i = iMin,iMax
443 CALL CHEAPAML_COARE3_FLUX(
444 I i, j, bi, bj, 0,
445 I theta(1-OLx,1-OLy,1,bi,bj), windSq,
446 O fsha0, flha0, evp_0, xolw0,
447 O ssqt0, q100_0, cdq_0,
448 O surfDrag(i,j,bi,bj),
449 O dumArg(1), dumArg(2), dumArg(3), dumArg(4),
450 I myIter, myThid )
451 Qnet(i,j,bi,bj) = (
452 & -solar(i,j,bi,bj)
453 & +xolw0 - cheapdlongwave(i,j,bi,bj)
454 & +fsha0
455 & +flha0
456 & )*maskC(i,j,1,bi,bj)
457 EmPmR(i,j,bi,bj) = evp_0
458 icFrac = iceFrac(i,j)
459 opFrac = 1. _d 0 - icFrac
460 C- Qsw (from FFIELDS.h) has opposite (and annoying) sign convention:
461 Qsw(i,j,bi,bj) = - ( icFrac*Qsw(i,j,bi,bj)
462 & + opFrac*solar(i,j,bi,bj) )
463 fsha(i,j) = icFrac*fsha(i,j) + opFrac*fsha0
464 flha(i,j) = icFrac*flha(i,j) + opFrac*flha0
465 evp(i,j) = icFrac*evp(i,j) + opFrac*evp_0
466 xolw(i,j) = icFrac*xolw(i,j) + opFrac*xolw0
467 ssqt(i,j) = icFrac*ssqt(i,j) + opFrac*ssqt0
468 q100(i,j) = icFrac*q100(i,j) + opFrac*q100_0
469 cdq(i,j) = icFrac*cdq(i,j) + opFrac*cdq_0
470 ENDDO
471 ENDDO
472 ELSE
473 DO j = jMin,jMax
474 DO i = iMin,iMax
475 IF (FluxFormula.EQ.'LANL') THEN
476 CALL cheapaml_LANL_flux(
477 I i, j, bi, bj,
478 O fsha(i,j), flha(i,j), evp(i,j),
479 O xolw(i,j), ssqt(i,j), q100(i,j) )
480 ELSEIF (FluxFormula.EQ.'COARE3') THEN
481 CALL CHEAPAML_COARE3_FLUX(
482 I i, j, bi, bj, 0,
483 I theta(1-OLx,1-OLy,1,bi,bj), windSq,
484 O fsha(i,j), flha(i,j), evp(i,j), xolw(i,j),
485 O ssqt(i,j), q100(i,j), cdq(i,j),
486 O surfDrag(i,j,bi,bj),
487 O dumArg(1), dumArg(2), dumArg(3), dumArg(4),
488 I myIter, myThid )
489 ENDIF
490 IF (useFreshWaterFlux) THEN
491 EmPmR(i,j,bi,bj) = evp(i,j)
492 ENDIF
493 ENDDO
494 ENDDO
495 ENDIF
496
497 DO j = jMin,jMax
498 DO i = iMin,iMax
499
500 C atmospheric upwelled long wave
501 ttt = Tair(i,j,bi,bj)-gamma_blk*(CheapHgrid(i,j,bi,bj)-zt)
502 c xalwu = stefan*(ttt+celsius2K)**4*0.5 _d 0
503 xalwu = stefan*(0.5*Tair(i,j,bi,bj)+0.5*ttt+celsius2K)**4
504 & *0.5 _d 0
505 C atmospheric downwelled long wave
506 xalwd = stefan*(Tair(i,j,bi,bj)+celsius2K)**4*0.5 _d 0
507 C total flux at upper atmospheric layer interface
508 xflu = ( -solar(i,j,bi,bj) + xalwu + flha(i,j)
509 & )*xef*maskC(i,j,1,bi,bj)
510 C lower flux calculation.
511 xfld = ( -solar(i,j,bi,bj) - xalwd + xolw(i,j)
512 & + fsha(i,j) + flha(i,j)
513 & )*xef*maskC(i,j,1,bi,bj)
514
515 IF (useDLongWave) THEN
516 xlwnet = xolw(i,j)-cheapdlongwave(i,j,bi,bj)
517 ELSE
518 C net long wave (see Josey et al. JGR 1997)
519 C coef lambda replaced by 0.5+lat/230
520 C convert spec humidity in water vapor pressure (mbar) using coef 1000/0.622=1607.7
521 xlwnet = 0.98 _d 0*stefan*(theta(i,j,1,bi,bj)+celsius2K)**4
522 & *(0.39 _d 0 - 0.05 _d 0*SQRT(qair(i,j,bi,bj)*1607.7 _d 0))
523 & *( oneRL - (halfRL+ABS(yG(i,j,bi,bj))/230. _d 0)
524 & *cheapclouds(i,j,bi,bj)**2 )
525 & + 4.0*0.98 _d 0*stefan*(theta(i,j,1,bi,bj)+celsius2K)**3
526 & *(theta(i,j,1,bi,bj)-Tair(i,j,bi,bj))
527
528 c xlwnet = xolw-stefan*(theta(i,j,1,bi,bj)+celsius2K)**4.
529 c & *(0.65+11.22*qair(i,j,bi,bj) + 0.25*cheapclouds(i,j,bi,bj)
530 c & -8.23*qair(i,j,bi,bj)*cheapclouds(i,j,bi,bj))
531 ENDIF
532 C clouds
533 c ttt2=Tair(i,j,bi,bj)-1.5*gamma_blk*CheapHgrid(i,j,bi,bj)
534 c Fclouds = stefan*ttt2**4*(0.4*cheapclouds(i,j,bi,bj)+1-0.4)/2
535 c ttt2=Tair(i,j,bi,bj)-3*gamma_blk*CheapHgrid(i,j,bi,bj)+celsius2K
536 c Fclouds = 0.3*stefan*ttt2**4 + 0.22*xolw*cheapclouds(i,j,bi,bj)
537 C add flux divergences into atmospheric temperature tendency
538 gTair(i,j,bi,bj)= (xfld-xflu)/CheapHgrid(i,j,bi,bj)
539 IF ( .NOT.useThSIce ) THEN
540 Qnet(i,j,bi,bj) = (
541 & -solar(i,j,bi,bj)
542 c & -xalwd
543 c & -Fclouds
544 c & +xolw
545 & +xlwnet
546 & +fsha(i,j)
547 & +flha(i,j)
548 & )*maskC(i,j,1,bi,bj)
549 Qsw(i,j,bi,bj) = -solar(i,j,bi,bj)
550 ENDIF
551
552 C need to precip?
553 IF (useFreshWaterFlux) THEN
554 q=q100(i,j)
555 C compute saturation specific humidity at atmospheric
556 C layer top
557 C first, what is the pressure there?
558 C ts is surface atmospheric temperature
559 ts=Tair(i,j,bi,bj)+gamma_blk*zt+celsius2K
560 pt=p0*(1-gamma_blk*CheapHgrid(i,j,bi,bj)/ts)
561 & **(gravity/gamma_blk/gasR)
562
563 C factor to compute rainfall from specific humidity
564 dm=100.*(p0-pt)*recip_gravity
565 C Large scale precip
566 precip = 0.
567 IF ( wWind(i,j,bi,bj).GT.0. .AND.
568 & q.GT.ssqt(i,j)*0.7 _d 0 ) THEN
569 precip = precip
570 & + ( (q-ssqt(i,j)*0.7 _d 0)*dm/cheap_pr2 )
571 & *(wWind(i,j,bi,bj)/0.75 _d -5)**2
572 ENDIF
573
574 C Convective precip
575 IF (q.GT.0.0214 _d 0 .AND. q.GT.ssqt(i,j)*0.9 _d 0) THEN
576 precip = precip + ((q-ssqt(i,j)*0.9 _d 0)*dm/cheap_pr1)
577 ENDIF
578
579 cheapPrecip(i,j,bi,bj) = precip*1200/CheapHgrid(i,j,bi,bj)
580 entrain = cdq(i,j)*q*0.25
581
582 c gqair(i,j,bi,bj)=(evp-precip-entrain)/CheapHgrid(i,j,bi,bj)
583 gqair(i,j,bi,bj) = (evp(i,j)-entrain)/CheapHgrid(i,j,bi,bj)
584 & /rhoa*maskC(i,j,1,bi,bj)
585 EmPmR(i,j,bi,bj) = ( EmPmR(i,j,bi,bj)
586 & -cheapPrecip(i,j,bi,bj)
587 & )*maskC(i,j,1,bi,bj)
588 ENDIF
589
590 ENDDO
591 ENDDO
592
593 C it is not necessary to use the Adams2d subroutine as
594 C the forcing is always computed at the current time step.
595 C note: full oceanic time step deltaT is used below
596 CALL CHEAPAML_TIMESTEP(
597 I bi, bj, iMin, iMax, jMin, jMax, deltaT,
598 I gTair,
599 U Tair,
600 I 0, myIter, myThid )
601 C do implicit time stepping over land
602 DO j=1-OLy,sNy+OLy
603 DO i=1-OLx,sNx+OLx
604 dtemp=tr(i,j,bi,bj)-Tair(i,j,bi,bj)
605 Tair(i,j,bi,bj)=Tair(i,j,bi,bj)+dtemp*xrelf(i,j,bi,bj)
606 ENDDO
607 ENDDO
608
609 C do water
610 IF (useFreshWaterFlux) THEN
611 CALL CHEAPAML_TIMESTEP(
612 I bi, bj, iMin, iMax, jMin, jMax, deltaT,
613 I gqair,
614 U qair,
615 I 0, myIter, myThid )
616 C do implicit time stepping over land and or buffer
617 DO j=1-OLy,sNy+OLy
618 DO i=1-OLx,sNx+OLx
619 dq=qr(i,j,bi,bj)-qair(i,j,bi,bj)
620 qair(i,j,bi,bj)=qair(i,j,bi,bj)+dq*xrelf(i,j,bi,bj)
621 IF (qair(i,j,bi,bj).LT.0.0) qair(i,j,bi,bj) = 0.0 _d 0
622 ENDDO
623 ENDDO
624 ENDIF
625
626 C do tracer
627 IF (useCheapTracer) THEN
628 C do implicit time stepping over land and or buffer
629 DO j=1-OLy,sNy+OLy
630 DO i=1-OLx,sNx+OLx
631 dtr=CheaptracerR(i,j,bi,bj)-Cheaptracer(i,j,bi,bj)
632 Cheaptracer(i,j,bi,bj) = Cheaptracer(i,j,bi,bj)
633 & + dtr*xrelf(i,j,bi,bj)
634 ENDDO
635 ENDDO
636 ENDIF
637
638 #ifdef ALLOW_DIAGNOSTICS
639 IF ( useDiagnostics ) THEN
640 CALL DIAGNOSTICS_FILL( fsha,'CH_SH ',0,1,2,bi,bj,myThid )
641 CALL DIAGNOSTICS_FILL( flha,'CH_LH ',0,1,2,bi,bj,myThid )
642 ENDIF
643 #endif /* ALLOW_DIAGNOSTICS */
644
645 C close bi,bj loops
646 ENDDO
647 ENDDO
648
649 C update edges
650 _EXCH_XY_RL(Tair,myThid)
651 _EXCH_XY_RS(Qnet,myThid)
652 IF (useFreshWaterFlux) THEN
653 _EXCH_XY_RL(qair,myThid)
654 _EXCH_XY_RS(EmPmR,myThid)
655 ENDIF
656 IF (useCheapTracer) THEN
657 _EXCH_XY_RL(Cheaptracer,myThid)
658 ENDIF
659 IF (.NOT.useStressOption) THEN
660 IF ( FluxFormula.EQ.'COARE3' ) THEN
661 _EXCH_XY_RL( surfDrag, myThid )
662 ELSE
663 CALL EXCH_UV_AGRID_3D_RL( ustress, vstress, .TRUE., 1, myThid )
664 ENDIF
665 ENDIF
666
667 C reset edges to open boundary profiles
668 c IF ( .NOT.(cheapamlXperiodic.AND.cheapamlYperiodic) ) THEN
669 IF ( notUsingXPeriodicity.OR.notUsingYPeriodicity ) THEN
670 xIsPeriodic = .NOT.notUsingXPeriodicity
671 yIsPeriodic = .NOT.notUsingYPeriodicity
672 DO bj=myByLo(myThid),myByHi(myThid)
673 DO bi=myBxLo(myThid),myBxHi(myThid)
674 CALL CHEAPAML_COPY_EDGES(
675 c I cheapamlXperiodic, cheapamlYperiodic,
676 I xIsPeriodic, yIsPeriodic,
677 I Tr(1-OLx,1-OLy,bi,bj),
678 U Tair(1-OLx,1-OLy,bi,bj),
679 I bi, bj, myIter, myThid )
680 IF (useFreshWaterFlux) THEN
681 CALL CHEAPAML_COPY_EDGES(
682 c I cheapamlXperiodic, cheapamlYperiodic,
683 I xIsPeriodic, yIsPeriodic,
684 I qr(1-OLx,1-OLy,bi,bj),
685 U qair(1-OLx,1-OLy,bi,bj),
686 I bi, bj, myIter, myThid )
687 ENDIF
688 IF (useCheapTracer) THEN
689 CALL CHEAPAML_COPY_EDGES(
690 c I cheapamlXperiodic, cheapamlYperiodic,
691 I xIsPeriodic, yIsPeriodic,
692 I CheaptracerR(1-OLx,1-OLy,bi,bj),
693 U Cheaptracer(1-OLx,1-OLy,bi,bj),
694 I bi, bj, myIter, myThid )
695 ENDIF
696 ENDDO
697 ENDDO
698 ENDIF
699
700 c CALL PLOT_FIELD_XYRS( gTair, 'S/R CHEAPAML gTair',1,myThid)
701 c CALL PLOT_FIELD_XYRS( Tair, 'S/R CHEAPAML Tair',1,myThid)
702 c CALL PLOT_FIELD_XYRS( Qnet, 'S/R CHEAPAML Qnet',1,myThid)
703
704 DO bj=myByLo(myThid),myByHi(myThid)
705 DO bi=myBxLo(myThid),myBxHi(myThid)
706
707 IF ( .NOT.useStressOption .AND. FluxFormula.EQ.'COARE3' ) THEN
708 DO j = 1-OLy+1,sNy+OLy
709 DO i = 1-OLx+1,sNx+OLx
710 fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
711 & *( surfDrag(i-1,j,bi,bj) + surfDrag(i,j,bi,bj) )
712 & *( uWind(i,j,bi,bj)-uVel(i,j,1,bi,bj) )
713 fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
714 & *( surfDrag(i,j-1,bi,bj) + surfDrag(i,j,bi,bj) )
715 & *( vWind(i,j,bi,bj)-vVel(i,j,1,bi,bj) )
716 ENDDO
717 ENDDO
718 #ifdef INCONSISTENT_WIND_LOCATION
719 DO j = 1-OLy,sNy+OLy
720 DO i = 1-OLx+1,sNx+OLx
721 fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
722 & *( surfDrag(i-1,j,bi,bj)
723 & *( uWind(i-1,j,bi,bj)-uVel(i-1,j,1,bi,bj) )
724 & + surfDrag(i,j,bi,bj)
725 & *( uWind(i,j,bi,bj) - uVel(i,j,1,bi,bj) ) )
726 ENDDO
727 ENDDO
728 DO j = 1-OLy+1,sNy+OLy
729 DO i = 1-OLx,sNx+OLx
730 fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
731 & *( surfDrag(i,j-1,bi,bj)
732 & *( vWind(i,j-1,bi,bj)-vVel(i,j-1,1,bi,bj) )
733 & + surfDrag(i,j,bi,bj)
734 & *( vWind(i,j,bi,bj) - vVel(i,j,1,bi,bj) ) )
735 ENDDO
736 ENDDO
737 #endif /* INCONSISTENT_WIND_LOCATION */
738 ELSE
739 Cswd move wind stresses to u and v points
740 DO j = 1-OLy,sNy+OLy
741 DO i = 1-OLx+1,sNx+OLx
742 fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
743 & *( ustress(i,j,bi,bj) + ustress(i-1,j,bi,bj) )
744 ENDDO
745 ENDDO
746 DO j = 1-OLy+1,sNy+OLy
747 DO i = 1-OLx,sNx+OLx
748 fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
749 & *( vstress(i,j,bi,bj) + vstress(i,j-1,bi,bj) )
750 ENDDO
751 ENDDO
752 ENDIF
753
754 C-- end bi,bj loops
755 ENDDO
756 ENDDO
757
758 #endif /* ALLOW_SEAGER */
759
760 #ifdef ALLOW_DIAGNOSTICS
761 IF ( useDiagnostics ) THEN
762 CALL DIAGNOSTICS_FILL(uWind, 'CH_Uwind',0,1,0,1,1,myThid)
763 CALL DIAGNOSTICS_FILL(vWind, 'CH_Vwind',0,1,0,1,1,myThid)
764 CALL DIAGNOSTICS_FILL_RS(Qnet,'CH_QNET ',0,1,0,1,1,myThid)
765 IF (useFreshWaterFlux) THEN
766 CALL DIAGNOSTICS_FILL_RS( EmPmR, 'CH_EmP ', 0,1,0,1,1,myThid)
767 CALL DIAGNOSTICS_FILL(cheapPrecip,'CH_Prec ',0,1,0,1,1,myThid)
768 ENDIF
769 IF (useCheapTracer) THEN
770 CALL DIAGNOSTICS_FILL(Cheaptracer,'CH_Trace',0,1,0,1,1,myThid)
771 ENDIF
772 ENDIF
773 #endif /* ALLOW_DIAGNOSTICS */
774
775 c DO bj=myByLo(myThid),myByHi(myThid)
776 c DO bi=myBxLo(myThid),myBxHi(myThid)
777 c DO j = 1-OLy,sNy+OLy
778 c DO i = 1-OLx+1,sNx+OLx
779 c fu(i,j,bi,bj) = 0.0
780 c fv(i,j,bi,bj) = 0.0
781 c Qnet(i,j,bi,bj) = 0.0
782 c EmPmR(i,j,bi,bj) = 0.0
783 c ENDDO
784 c ENDDO
785 c ENDDO
786 c ENDDO
787
788 RETURN
789 END
  ViewVC Help
Powered by ViewVC 1.1.22