/[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.1 - (hide annotations) (download)
Wed May 22 19:39:51 2013 UTC (12 years, 1 month ago) by jmc
Branch: MAIN
unfinished set-up to test pkg/cheapaml with thsice

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

  ViewVC Help
Powered by ViewVC 1.1.22