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

Annotation 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.4 - (hide annotations) (download)
Mon Jun 17 13:48:55 2013 UTC (12 years ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +1 -1 lines
FILE REMOVED
source files now part of standard pkg/cheapaml code

1 jmc 1.4 C $Header: /u/gcmpack/MITgcm_contrib/verification_other/offline_cheapaml/code/cheapaml.F,v 1.3 2013/06/11 01:57:46 jmc Exp $
2 jmc 1.1 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 jmc 1.2 C uWind, vWind - mean wind speed (m/s)
28 jmc 1.1 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 jmc 1.2 INTEGER myThid
51 jmc 1.1
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 jmc 1.2 _RL windSq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77 jmc 1.1 _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 jmc 1.2 C surfDrag :: surface drag coeff (for wind stress)
85     _RL surfDrag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
86 jmc 1.1 _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 jmc 1.2 uTrans(i,j)=uWind(i,j,bi,bj)
259     vTrans(i,j)=vWind(i,j,bi,bj)
260 jmc 1.1 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 jmc 1.2 I uWind, vWind,
266 jmc 1.1 I cheapaml_kdiff, Tair,
267     I deltaTtracer, zu, useFluxLimit,
268     I cheapamlXperiodic, cheapamlYperiodic,
269 jmc 1.2 O wWind,
270 jmc 1.1 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 jmc 1.2 _EXCH_XY_RL(Tair,myThid)
287 jmc 1.1
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 jmc 1.2 uTrans(i,j)=uWind(i,j,bi,bj)
296     vTrans(i,j)=vWind(i,j,bi,bj)
297 jmc 1.1 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 jmc 1.2 I uWind, vWind,
303 jmc 1.1 I cheapaml_kdiff, qair,
304     I deltaTtracer, zu, useFluxLimit,
305     I cheapamlXperiodic, cheapamlYperiodic,
306 jmc 1.2 O wWind,
307 jmc 1.1 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 jmc 1.2 _EXCH_XY_RL(qair,myThid)
324 jmc 1.1 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 jmc 1.2 uTrans(i,j)=uWind(i,j,bi,bj)
334     vTrans(i,j)=vWind(i,j,bi,bj)
335 jmc 1.1 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 jmc 1.2 I uWind, vWind,
341 jmc 1.1 I cheapaml_kdiff, Cheaptracer,
342     I deltaTtracer, zu, useFluxLimit,
343     I cheapamlXperiodic, cheapamlYperiodic,
344 jmc 1.2 O wWind,
345 jmc 1.1 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 jmc 1.2 _EXCH_XY_RL(Cheaptracer,myThid)
362 jmc 1.1 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 jmc 1.2 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 jmc 1.1 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 jmc 1.2 I uWind(1-OLx,1-OLy,bi,bj),
437     I vWind(1-OLx,1-OLy,bi,bj), lath,
438 jmc 1.1 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 jmc 1.2 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 jmc 1.1 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 jmc 1.2 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 jmc 1.1 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 jmc 1.2 ttt = Tair(i,j,bi,bj)-gamma_blk*(CheapHgrid(i,j,bi,bj)-zt)
502 jmc 1.1 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 jmc 1.2 & *0.5 _d 0
505 jmc 1.1 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 jmc 1.2 c xlwnet = xolw-stefan*(theta(i,j,1,bi,bj)+celsius2K)**4.
529 jmc 1.1 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 jmc 1.2 c ttt2=Tair(i,j,bi,bj)-1.5*gamma_blk*CheapHgrid(i,j,bi,bj)
534 jmc 1.1 c Fclouds = stefan*ttt2**4*(0.4*cheapclouds(i,j,bi,bj)+1-0.4)/2
535 jmc 1.2 c ttt2=Tair(i,j,bi,bj)-3*gamma_blk*CheapHgrid(i,j,bi,bj)+celsius2K
536 jmc 1.1 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 jmc 1.2 gTair(i,j,bi,bj)= (xfld-xflu)/CheapHgrid(i,j,bi,bj)
539 jmc 1.1 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 jmc 1.2 pt=p0*(1-gamma_blk*CheapHgrid(i,j,bi,bj)/ts)
561 jmc 1.1 & **(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 jmc 1.2 IF ( wWind(i,j,bi,bj).GT.0. .AND.
568 jmc 1.1 & 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 jmc 1.2 & *(wWind(i,j,bi,bj)/0.75 _d -5)**2
572 jmc 1.1 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 jmc 1.2 cheapPrecip(i,j,bi,bj) = precip*1200/CheapHgrid(i,j,bi,bj)
580 jmc 1.1 entrain = cdq(i,j)*q*0.25
581    
582 jmc 1.2 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 jmc 1.1 & /rhoa*maskC(i,j,1,bi,bj)
585     EmPmR(i,j,bi,bj) = ( EmPmR(i,j,bi,bj)
586 jmc 1.2 & -cheapPrecip(i,j,bi,bj)
587 jmc 1.1 & )*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 jmc 1.3 CALL DIAGNOSTICS_FILL( q100,'CH_q100 ',0,1,2,bi,bj,myThid )
643     CALL DIAGNOSTICS_FILL( ssqt,'CH_ssqt ',0,1,2,bi,bj,myThid )
644 jmc 1.1 ENDIF
645     #endif /* ALLOW_DIAGNOSTICS */
646    
647     C close bi,bj loops
648     ENDDO
649     ENDDO
650    
651     C update edges
652 jmc 1.2 _EXCH_XY_RL(Tair,myThid)
653     _EXCH_XY_RS(Qnet,myThid)
654 jmc 1.1 IF (useFreshWaterFlux) THEN
655 jmc 1.2 _EXCH_XY_RL(qair,myThid)
656     _EXCH_XY_RS(EmPmR,myThid)
657 jmc 1.1 ENDIF
658     IF (useCheapTracer) THEN
659 jmc 1.2 _EXCH_XY_RL(Cheaptracer,myThid)
660 jmc 1.1 ENDIF
661     IF (.NOT.useStressOption) THEN
662 jmc 1.2 IF ( FluxFormula.EQ.'COARE3' ) THEN
663     _EXCH_XY_RL( surfDrag, myThid )
664     ELSE
665     CALL EXCH_UV_AGRID_3D_RL( ustress, vstress, .TRUE., 1, myThid )
666     ENDIF
667 jmc 1.1 ENDIF
668    
669     C reset edges to open boundary profiles
670     c IF ( .NOT.(cheapamlXperiodic.AND.cheapamlYperiodic) ) THEN
671     IF ( notUsingXPeriodicity.OR.notUsingYPeriodicity ) THEN
672     xIsPeriodic = .NOT.notUsingXPeriodicity
673     yIsPeriodic = .NOT.notUsingYPeriodicity
674     DO bj=myByLo(myThid),myByHi(myThid)
675     DO bi=myBxLo(myThid),myBxHi(myThid)
676     CALL CHEAPAML_COPY_EDGES(
677     c I cheapamlXperiodic, cheapamlYperiodic,
678     I xIsPeriodic, yIsPeriodic,
679     I Tr(1-OLx,1-OLy,bi,bj),
680     U Tair(1-OLx,1-OLy,bi,bj),
681     I bi, bj, myIter, myThid )
682     IF (useFreshWaterFlux) THEN
683     CALL CHEAPAML_COPY_EDGES(
684     c I cheapamlXperiodic, cheapamlYperiodic,
685     I xIsPeriodic, yIsPeriodic,
686     I qr(1-OLx,1-OLy,bi,bj),
687     U qair(1-OLx,1-OLy,bi,bj),
688     I bi, bj, myIter, myThid )
689     ENDIF
690     IF (useCheapTracer) THEN
691     CALL CHEAPAML_COPY_EDGES(
692     c I cheapamlXperiodic, cheapamlYperiodic,
693     I xIsPeriodic, yIsPeriodic,
694     I CheaptracerR(1-OLx,1-OLy,bi,bj),
695     U Cheaptracer(1-OLx,1-OLy,bi,bj),
696     I bi, bj, myIter, myThid )
697     ENDIF
698     ENDDO
699     ENDDO
700     ENDIF
701    
702     c CALL PLOT_FIELD_XYRS( gTair, 'S/R CHEAPAML gTair',1,myThid)
703     c CALL PLOT_FIELD_XYRS( Tair, 'S/R CHEAPAML Tair',1,myThid)
704     c CALL PLOT_FIELD_XYRS( Qnet, 'S/R CHEAPAML Qnet',1,myThid)
705    
706     DO bj=myByLo(myThid),myByHi(myThid)
707     DO bi=myBxLo(myThid),myBxHi(myThid)
708 jmc 1.2
709     IF ( .NOT.useStressOption .AND. FluxFormula.EQ.'COARE3' ) THEN
710     DO j = 1-OLy+1,sNy+OLy
711     DO i = 1-OLx+1,sNx+OLx
712     fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
713     & *( surfDrag(i-1,j,bi,bj) + surfDrag(i,j,bi,bj) )
714     & *( uWind(i,j,bi,bj)-uVel(i,j,1,bi,bj) )
715     fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
716     & *( surfDrag(i,j-1,bi,bj) + surfDrag(i,j,bi,bj) )
717     & *( vWind(i,j,bi,bj)-vVel(i,j,1,bi,bj) )
718     ENDDO
719     ENDDO
720     #ifdef INCONSISTENT_WIND_LOCATION
721     DO j = 1-OLy,sNy+OLy
722     DO i = 1-OLx+1,sNx+OLx
723     fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
724     & *( surfDrag(i-1,j,bi,bj)
725     & *( uWind(i-1,j,bi,bj)-uVel(i-1,j,1,bi,bj) )
726     & + surfDrag(i,j,bi,bj)
727     & *( uWind(i,j,bi,bj) - uVel(i,j,1,bi,bj) ) )
728     ENDDO
729     ENDDO
730     DO j = 1-OLy+1,sNy+OLy
731     DO i = 1-OLx,sNx+OLx
732     fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
733     & *( surfDrag(i,j-1,bi,bj)
734     & *( vWind(i,j-1,bi,bj)-vVel(i,j-1,1,bi,bj) )
735     & + surfDrag(i,j,bi,bj)
736     & *( vWind(i,j,bi,bj) - vVel(i,j,1,bi,bj) ) )
737     ENDDO
738     ENDDO
739     #endif /* INCONSISTENT_WIND_LOCATION */
740     ELSE
741 jmc 1.1 Cswd move wind stresses to u and v points
742     DO j = 1-OLy,sNy+OLy
743     DO i = 1-OLx+1,sNx+OLx
744 jmc 1.2 fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
745     & *( ustress(i,j,bi,bj) + ustress(i-1,j,bi,bj) )
746 jmc 1.1 ENDDO
747     ENDDO
748     DO j = 1-OLy+1,sNy+OLy
749     DO i = 1-OLx,sNx+OLx
750 jmc 1.2 fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
751     & *( vstress(i,j,bi,bj) + vstress(i,j-1,bi,bj) )
752 jmc 1.1 ENDDO
753     ENDDO
754 jmc 1.2 ENDIF
755 jmc 1.1
756     C-- end bi,bj loops
757     ENDDO
758     ENDDO
759    
760     #endif /* ALLOW_SEAGER */
761    
762     #ifdef ALLOW_DIAGNOSTICS
763     IF ( useDiagnostics ) THEN
764 jmc 1.2 CALL DIAGNOSTICS_FILL(uWind, 'CH_Uwind',0,1,0,1,1,myThid)
765     CALL DIAGNOSTICS_FILL(vWind, 'CH_Vwind',0,1,0,1,1,myThid)
766 jmc 1.1 CALL DIAGNOSTICS_FILL_RS(Qnet,'CH_QNET ',0,1,0,1,1,myThid)
767 jmc 1.2 IF (useFreshWaterFlux) THEN
768     CALL DIAGNOSTICS_FILL_RS( EmPmR, 'CH_EmP ', 0,1,0,1,1,myThid)
769     CALL DIAGNOSTICS_FILL(cheapPrecip,'CH_Prec ',0,1,0,1,1,myThid)
770     ENDIF
771 jmc 1.1 IF (useCheapTracer) THEN
772     CALL DIAGNOSTICS_FILL(Cheaptracer,'CH_Trace',0,1,0,1,1,myThid)
773     ENDIF
774     ENDIF
775     #endif /* ALLOW_DIAGNOSTICS */
776    
777     c DO bj=myByLo(myThid),myByHi(myThid)
778     c DO bi=myBxLo(myThid),myBxHi(myThid)
779     c DO j = 1-OLy,sNy+OLy
780     c DO i = 1-OLx+1,sNx+OLx
781     c fu(i,j,bi,bj) = 0.0
782     c fv(i,j,bi,bj) = 0.0
783     c Qnet(i,j,bi,bj) = 0.0
784     c EmPmR(i,j,bi,bj) = 0.0
785     c ENDDO
786     c ENDDO
787     c ENDDO
788     c ENDDO
789    
790     RETURN
791     END

  ViewVC Help
Powered by ViewVC 1.1.22