ROMS
Loading...
Searching...
No Matches
ad_bulk_flux.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE ad_bulk_flux_mod
3#if defined ADJOINT && defined BULK_FLUXES_NOT_YET
4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11! !
12! This routine computes the bulk parameterization of surface wind !
13! stress and surface net heat fluxes. !
14! !
15! References: !
16! !
17! Fairall, C.W., E.F. Bradley, D.P. Rogers, J.B. Edson and G.S. !
18! Young, 1996: Bulk parameterization of air-sea fluxes for !
19! tropical ocean-global atmosphere Coupled-Ocean Atmosphere !
20! Response Experiment, JGR, 101, 3747-3764. !
21! !
22! Fairall, C.W., E.F. Bradley, J.S. Godfrey, G.A. Wick, J.B. !
23! Edson, and G.S. Young, 1996: Cool-skin and warm-layer !
24! effects on sea surface temperature, JGR, 101, 1295-1308. !
25! !
26! Liu, W.T., K.B. Katsaros, and J.A. Businger, 1979: Bulk !
27! parameterization of the air-sea exchange of heat and !
28! water vapor including the molecular constraints at !
29! the interface, J. Atmos. Sci, 36, 1722-1735. !
30! !
31! Adapted from COARE code written originally by David Rutgers and !
32! Frank Bradley. !
33! !
34! EMINUSP option for equivalent salt fluxes added by Paul Goodman !
35! (10/2004). !
36! !
37! Modified by Kate Hedstrom for COARE version 3.0 (03/2005). !
38! Modified by Jim Edson to correct specific hunidities. !
39! !
40! Reference: !
41! !
42! Fairall et al., 2003: J. Climate, 16, 571-591. !
43! !
44!=======================================================================
45!
46 implicit none
47!
48 PRIVATE
49 PUBLIC :: ad_bulk_flux
50!
51 CONTAINS
52!
53!***********************************************************************
54 SUBROUTINE ad_bulk_flux (ng, tile)
55!***********************************************************************
56!
57 USE mod_param
58 USE mod_forces
59 USE mod_grid
60 USE mod_mixing
61 USE mod_ocean
62 USE mod_stepping
63!
64! Imported variable declarations.
65!
66 integer, intent(in) :: ng, tile
67!
68! Local variable declarations.
69!
70 character (len=*), parameter :: myfile = &
71 & __FILE__
72!
73# include "tile.h"
74!
75# ifdef PROFILE
76 CALL wclock_on (ng, iadm, 17, __line__, myfile)
77# endif
78 CALL ad_bulk_flux_tile (ng, tile, &
79 & lbi, ubi, lbj, ubj, &
80 & imins, imaxs jmins, jmaxs, &
81 & nrhs(ng), &
82# ifdef MASKING
83 & grid(ng) % rmask, &
84 & grid(ng) % umask, &
85 & grid(ng) % vmask, &
86# endif
87 & mixing(ng) % alpha, &
88 & mixing(ng) % ad_alpha, &
89 & mixing(ng) % beta, &
90 & mixing(ng) % ad_beta, &
91 & ocean(ng) % rho, &
92 & ocean(ng) % ad_rho, &
93 & ocean(ng) % t, &
94 & ocean(ng) % ad_t, &
95 & forces(ng) % Hair, &
96 & forces(ng) % Pair, &
97 & forces(ng) % Tair, &
98 & forces(ng) % Uwind, &
99 & forces(ng) % Vwind, &
100# ifdef CLOUDS
101 & forces(ng) % cloud, &
102# endif
103# ifdef BBL_MODEL_NOT_YET
104 & forces(ng) % Awave, &
105 & forces(ng) % Pwave, &
106# endif
107 & forces(ng) % rain, &
108 & forces(ng) % lhflx, &
109 & forces(ng) % ad_lhflx, &
110 & forces(ng) % lrflx, &
111 & forces(ng) % ad_lrflx, &
112 & forces(ng) % shflx, &
113 & forces(ng) % ad_shflx, &
114 & forces(ng) % srflx, &
115 & forces(ng) % stflx, &
116 & forces(ng) % ad_stflx, &
117# ifdef EMINUSP
118 & forces(ng) % evap, &
119 & forces(ng) % ad_evap, &
120# endif
121 & forces(ng) % ad_sustr, &
122 & forces(ng) % ad_svstr)
123
124# ifdef PROFILE
125 CALL wclock_off (ng, iadm, 17, __line__, myfile)
126# endif
127!
128 RETURN
129 END SUBROUTINE ad_bulk_flux
130!
131!***********************************************************************
132 SUBROUTINE ad_bulk_flux_tile (ng, tile, &
133 & LBi, UBi, LBj, UBj, &
134 & IminS, ImaxS, JminS, JmaxS &
135 & nrhs, &
136# ifdef MASKING
137 & rmask, umask, vmask, &
138# endif
139 & alpha, ad_alpha, &
140 & beta, ad_beta, &
141 & rho, ad_rho, t, ad_t, &
142 & Hair, Pair, Tair, &
143 & Uwind, Vwind, &
144# ifdef CLOUDS
145 & cloud, &
146# endif
147# ifdef BBL_MODEL_NOT_YET
148 & Awave, Pwave, &
149# endif
150 & rain, &
151 & lhflx, ad_lhflx, &
152 & lrflx, ad_lrflx, &
153 & shflx, ad_shflx, &
154 & srflx, &
155 & stflx, ad_stflx, &
156# ifdef EMINUSP
157 & evap, ad_evap, &
158# endif
159 & ad_sustr, ad_svstr)
160!***********************************************************************
161!
162 USE mod_param
163 USE mod_scalars
164!
165 USE bulk_flux_mod, ONLY : bulk_psiu, bulk_psit
166 USE exchange_2d_mod
167# ifdef DISTRIBUTE
168 USE mp_exchange_mod, ONLY : ad_mp_exchange2d
169# endif
170!
171! Imported variable declarations.
172!
173 integer, intent(in) :: ng, tile
174 integer, intent(in) :: LBi, UBi, LBj, UBj
175 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
176 integer, intent(in) :: nrhs
177!
178# ifdef ASSUMED_SHAPE
179# ifdef MASKING
180 real(r8), intent(in) :: rmask(LBi:,LBj:)
181 real(r8), intent(in) :: umask(LBi:,LBj:)
182 real(r8), intent(in) :: vmask(LBi:,LBj:)
183# endif
184 real(r8), intent(in) :: alpha(LBi:,LBj:)
185 real(r8), intent(in) :: beta(LBi:,LBj:)
186 real(r8), intent(in) :: rho(LBi:,LBj:,:)
187 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
188 real(r8), intent(in) :: Hair(LBi:,LBj:)
189 real(r8), intent(in) :: Pair(LBi:,LBj:)
190 real(r8), intent(in) :: Tair(LBi:,LBj:)
191 real(r8), intent(in) :: Uwind(LBi:,LBj:)
192 real(r8), intent(in) :: Vwind(LBi:,LBj:)
193 real(r8), intent(inout) :: ad_alpha(LBi:,LBj:)
194 real(r8), intent(inout) :: ad_beta(LBi:,LBj:)
195 real(r8), intent(inout) :: ad_rho(LBi:,LBj:,:)
196 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
197# ifdef CLOUDS
198 real(r8), intent(in) :: cloud(LBi:,LBj:)
199# endif
200# ifdef BBL_MODEL_NOT_YET
201 real(r8), intent(in) :: Awave(LBi:,LBj:)
202 real(r8), intent(in) :: Pwave(LBi:,LBj:)
203# endif
204 real(r8), intent(in) :: rain(LBi:,LBj:)
205
206 real(r8), intent(inout) :: lhflx(LBi:,LBj:)
207 real(r8), intent(inout) :: lrflx(LBi:,LBj:)
208 real(r8), intent(inout) :: shflx(LBi:,LBj:)
209 real(r8), intent(inout) :: srflx(LBi:,LBj:)
210 real(r8), intent(inout) :: stflx(LBi:,LBj:,:)
211 real(r8), intent(inout) :: ad_lhflx(LBi:,LBj:)
212 real(r8), intent(inout) :: ad_lrflx(LBi:,LBj:)
213 real(r8), intent(inout) :: ad_shflx(LBi:,LBj:)
214 real(r8), intent(inout) :: ad_stflx(LBi:,LBj:,:)
215
216# ifdef EMINUSP
217 real(r8), intent(in) :: evap(LBi:,LBj:)
218 real(r8), intent(inout) :: ad_evap(LBi:,LBj:)
219# endif
220 real(r8), intent(inout) :: ad_sustr(LBi:,LBj:)
221 real(r8), intent(inout) :: ad_svstr(LBi:,LBj:)
222# else
223# ifdef MASKING
224 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
225 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
226 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
227# endif
228 real(r8), intent(in) :: alpha(LBi:UBi,LBj:UBj)
229 real(r8), intent(in) :: beta(LBi:UBi,LBj:UBj)
230 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
231 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
232 real(r8), intent(in) :: Hair(LBi:UBi,LBj:UBj)
233 real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
234 real(r8), intent(in) :: Tair(LBi:UBi,LBj:UBj)
235 real(r8), intent(in) :: Uwind(LBi:UBi,LBj:UBj)
236 real(r8), intent(in) :: Vwind(LBi:UBi,LBj:UBj)
237 real(r8), intent(inout) :: ad_alpha(LBi:UBi,LBj:UBj)
238 real(r8), intent(inout) :: ad_beta(LBi:UBi,LBj:UBj)
239 real(r8), intent(inout) :: ad_rho(LBi:UBi,LBj:UBj,N(ng))
240 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
241# ifdef CLOUDS
242 real(r8), intent(in) :: cloud(LBi:UBi,LBj:UBj)
243# endif
244# ifdef BBL_MODEL_NOT_YET
245 real(r8), intent(in) :: Awave(LBi:UBi,LBj:UBj)
246 real(r8), intent(in) :: Pwave(LBi:UBi,LBj:UBj)
247# endif
248 real(r8), intent(in) :: rain(LBi:UBi,LBj:UBj)
249
250 real(r8), intent(inout) :: lhflx(LBi:UBi,LBj:UBj)
251 real(r8), intent(inout) :: lrflx(LBi:UBi,LBj:UBj)
252 real(r8), intent(inout) :: shflx(LBi:UBi,LBj:UBj)
253 real(r8), intent(inout) :: srflx(LBi:UBi,LBj:UBj)
254 real(r8), intent(inout) :: stflx(LBi:UBi,LBj:UBj,NT(ng))
255 real(r8), intent(inout) :: ad_lhflx(LBi:UBi,LBj:UBj)
256 real(r8), intent(inout) :: ad_lrflx(LBi:UBi,LBj:UBj)
257 real(r8), intent(inout) :: ad_shflx(LBi:UBi,LBj:UBj)
258 real(r8), intent(inout) :: ad_stflx(LBi:UBi,LBj:UBj,NT(ng))
259
260# ifdef EMINUSP
261 real(r8), intent(in) :: evap(LBi:UBi,LBj:UBj)
262 real(r8), intent(inout) :: ad_evap(LBi:UBi,LBj:UBj)
263# endif
264 real(r8), intent(inout) :: ad_sustr(LBi:UBi,LBj:UBj)
265 real(r8), intent(inout) :: ad_svstr(LBi:UBi,LBj:UBj)
266# endif
267!
268! Local variable declarations.
269!
270 integer :: Iter, IterA, i, j, k
271 integer, parameter :: IterMax = 3
272
273 real(r8), parameter :: eps = 1.0e-20_r8
274 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
275
276 real(r8) :: Bf, Cd, Hl, Hlw, Hscale, Hs, Hsr, IER
277 real(r8) :: ad_Bf, ad_Cd, ad_Hl, ad_Hlw, ad_Hsr, ad_Hs
278 real(r8) :: PairM, RH, Taur
279 real(r8) :: Wspeed, ZQoL, ZToL
280
281 real(r8) :: cff, cff1, cff2, diffh, diffw, oL, upvel
282 real(r8) :: twopi_inv, wet_bulb
283 real(r8) :: ad_wet_bulb, ad_Wspeed
284 real(r8) :: fac, ad_fac, fac1, ad_fac1, fac2, ad_fac2
285 real(r8) :: fac3, ad_fac3, ad_cff, ad_upvel
286 real(r8) :: adfac, adfac1, adfac2
287# ifdef LONGWAVE
288 real(r8) :: e_sat, vap_p
289# endif
290# ifdef COOL_SKIN
291 real(r8) :: Clam, Fc, Hcool, Hsb, Hlb, Qbouy, Qcool, lambd
292 real(r8) :: ad_Clam, ad_Hcool, ad_Hsb, ad_Hlb
293 real(r8) :: ad_Qbouy, ad_Qcool, ad_lambd
294# endif
295
296 real(r8), dimension(IminS:ImaxS) :: CC
297 real(r8), dimension(IminS:ImaxS) :: Cd10
298 real(r8), dimension(IminS:ImaxS) :: Ch10
299 real(r8), dimension(IminS:ImaxS) :: Ct10
300 real(r8), dimension(IminS:ImaxS) :: charn
301 real(r8), dimension(IminS:ImaxS) :: Ct
302 real(r8), dimension(IminS:ImaxS) :: Cwave
303 real(r8), dimension(IminS:ImaxS) :: Cwet
304 real(r8), dimension(IminS:ImaxS) :: delQ
305 real(r8), dimension(IminS:ImaxS) :: delQc
306 real(r8), dimension(IminS:ImaxS) :: delT
307 real(r8), dimension(IminS:ImaxS) :: delTc
308 real(r8), dimension(IminS:ImaxS) :: delW
309 real(r8), dimension(IminS:ImaxS) :: delW1
310 real(r8), dimension(IminS:ImaxS) :: L
311 real(r8), dimension(IminS:ImaxS) :: L10
312 real(r8), dimension(IminS:ImaxS) :: Lwave
313 real(r8), dimension(IminS:ImaxS) :: Q
314 real(r8), dimension(IminS:ImaxS) :: Qair
315 real(r8), dimension(IminS:ImaxS) :: Qpsi
316 real(r8), dimension(IminS:ImaxS) :: Qsea
317 real(r8), dimension(IminS:ImaxS) :: Qstar
318 real(r8), dimension(IminS:ImaxS) :: rhoAir
319 real(r8), dimension(IminS:ImaxS) :: rhoSea
320 real(r8), dimension(IminS:ImaxS) :: Ri
321 real(r8), dimension(IminS:ImaxS) :: Ribcu
322 real(r8), dimension(IminS:ImaxS) :: Rr
323 real(r8), dimension(IminS:ImaxS) :: Scff
324 real(r8), dimension(IminS:ImaxS) :: TairC
325 real(r8), dimension(IminS:ImaxS) :: TairK
326 real(r8), dimension(IminS:ImaxS) :: Tcff
327 real(r8), dimension(IminS:ImaxS) :: Tpsi
328 real(r8), dimension(IminS:ImaxS) :: TseaC
329 real(r8), dimension(IminS:ImaxS) :: TseaK
330 real(r8), dimension(IminS:ImaxS) :: Tstar
331 real(r8), dimension(IminS:ImaxS) :: u10
332 real(r8), dimension(IminS:ImaxS) :: VisAir
333 real(r8), dimension(IminS:ImaxS) :: Wgus
334 real(r8), dimension(IminS:ImaxS) :: Wmag
335 real(r8), dimension(IminS:ImaxS) :: Wpsi
336 real(r8), dimension(IminS:ImaxS) :: Wstar
337 real(r8), dimension(IminS:ImaxS) :: Zetu
338 real(r8), dimension(IminS:ImaxS) :: Zo10
339 real(r8), dimension(IminS:ImaxS) :: ZoT10
340 real(r8), dimension(IminS:ImaxS) :: ZoL
341 real(r8), dimension(IminS:ImaxS) :: ZoQ
342 real(r8), dimension(IminS:ImaxS) :: ZoT
343 real(r8), dimension(IminS:ImaxS) :: ZoW
344 real(r8), dimension(IminS:ImaxS) :: ZWoL
345
346 real(r8), dimension(IminS:ImaxS) :: Wstar1
347 real(r8), dimension(IminS:ImaxS) :: Tstar1
348 real(r8), dimension(IminS:ImaxS) :: Qstar1
349 real(r8), dimension(IminS:ImaxS) :: delTc1
350
351 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hlv
352 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LHeat
353 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LRad
354 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: SHeat
355 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: SRad
356 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taux
357 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauy
358
359 real(r8), dimension(IminS:ImaxS) :: ad_CC
360 real(r8), dimension(IminS:ImaxS) :: ad_charn
361 real(r8), dimension(IminS:ImaxS) :: ad_Ct
362 real(r8), dimension(IminS:ImaxS) :: ad_Cd10
363 real(r8), dimension(IminS:ImaxS) :: ad_Ct10
364 real(r8), dimension(IminS:ImaxS) :: ad_Cwet
365 real(r8), dimension(IminS:ImaxS) :: ad_delQ
366 real(r8), dimension(IminS:ImaxS) :: ad_delQc
367 real(r8), dimension(IminS:ImaxS) :: ad_delT
368 real(r8), dimension(IminS:ImaxS) :: ad_delTc
369 real(r8), dimension(IminS:ImaxS) :: ad_delW
370 real(r8), dimension(IminS:ImaxS) :: ad_L
371 real(r8), dimension(IminS:ImaxS) :: ad_L10
372 real(r8), dimension(IminS:ImaxS) :: ad_Q
373 real(r8), dimension(IminS:ImaxS) :: ad_Qpsi
374 real(r8), dimension(IminS:ImaxS) :: ad_Qsea
375 real(r8), dimension(IminS:ImaxS) :: ad_Qstar
376 real(r8), dimension(IminS:ImaxS) :: ad_rhoSea
377 real(r8), dimension(IminS:ImaxS) :: ad_Ri
378 real(r8), dimension(IminS:ImaxS) :: ad_Rr
379 real(r8), dimension(IminS:ImaxS) :: ad_Scff
380 real(r8), dimension(IminS:ImaxS) :: ad_Tcff
381 real(r8), dimension(IminS:ImaxS) :: ad_Tpsi
382 real(r8), dimension(IminS:ImaxS) :: ad_TseaC
383 real(r8), dimension(IminS:ImaxS) :: ad_TseaK
384 real(r8), dimension(IminS:ImaxS) :: ad_Tstar
385 real(r8), dimension(IminS:ImaxS) :: ad_u10
386 real(r8), dimension(IminS:ImaxS) :: ad_Wgus
387 real(r8), dimension(IminS:ImaxS) :: ad_Wpsi
388 real(r8), dimension(IminS:ImaxS) :: ad_Wstar
389 real(r8), dimension(IminS:ImaxS) :: ad_Zetu
390 real(r8), dimension(IminS:ImaxS) :: ad_Zo10
391 real(r8), dimension(IminS:ImaxS) :: ad_ZoT10
392 real(r8), dimension(IminS:ImaxS) :: ad_ZoL
393 real(r8), dimension(IminS:ImaxS) :: ad_ZoQ
394 real(r8), dimension(IminS:ImaxS) :: ad_ZoT
395 real(r8), dimension(IminS:ImaxS) :: ad_ZoW
396 real(r8), dimension(IminS:ImaxS) :: ad_ZWoL
397
398 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_Hlv
399 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_LHeat
400 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_LRad
401 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_SHeat
402 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_Taux
403 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_Tauy
404
405# include "set_bounds.h"
406!
407!------------------------------------------------------------------------
408! Initialize adjoint private variables.
409!------------------------------------------------------------------------
410!
411 ad_wet_bulb=0.0_r8
412 ad_wspeed=0.0_r8
413 ad_bf=0.0_r8
414 ad_cd=0.0_r8
415 ad_hl=0.0_r8
416 ad_hlw=0.0_r8
417 ad_hsr=0.0_r8
418 ad_hs=0.0_r8
419 ad_fac=0.0_r8
420 ad_fac1=0.0_r8
421 ad_fac2=0.0_r8
422 ad_fac3=0.0_r8
423 ad_cff=0.0_r8
424 ad_upvel=0.0_r8
425# ifdef COOL_SKIN
426 ad_clam=0.0_r8
427 ad_hcool=0.0_r8
428 ad_hsb=0.0_r8
429 ad_hlb=0.0_r8
430 ad_qbouy=0.0_r8
431 ad_qcool=0.0_r8
432 ad_lambd=0.0_r8
433# endif
434 DO i=imins,imaxs
435 ad_cc(i)=0.0_r8
436 ad_charn(i)=0.0_r8
437 ad_ct(i)=0.0_r8
438 ad_cd10(i)=0.0_r8
439 ad_ct10(i)=0.0_r8
440 ad_cwet(i)=0.0_r8
441 ad_delq(i)=0.0_r8
442 ad_delqc(i)=0.0_r8
443 ad_delt(i)=0.0_r8
444 ad_deltc(i)=0.0_r8
445 ad_delw(i)=0.0_r8
446 ad_l(i)=0.0_r8
447 ad_l10(i)=0.0_r8
448 ad_q(i)=0.0_r8
449 ad_qpsi(i)=0.0_r8
450 ad_qsea(i)=0.0_r8
451 ad_qstar(i)=0.0_r8
452 ad_rhosea(i)=0.0_r8
453 ad_ri(i)=0.0_r8
454 ad_rr(i)=0.0_r8
455 ad_scff(i)=0.0_r8
456 ad_tcff(i)=0.0_r8
457 ad_tpsi(i)=0.0_r8
458 ad_tseac(i)=0.0_r8
459 ad_tseak(i)=0.0_r8
460 ad_tstar(i)=0.0_r8
461 ad_u10(i)=0.0_r8
462 ad_wgus(i)=0.0_r8
463 ad_wpsi(i)=0.0_r8
464 ad_wstar(i)=0.0_r8
465 ad_zetu(i)=0.0_r8
466 ad_zo10(i)=0.0_r8
467 ad_zot10(i)=0.0_r8
468 ad_zol(i)=0.0_r8
469 ad_zoq(i)=0.0_r8
470 ad_zot(i)=0.0_r8
471 ad_zow(i)=0.0_r8
472 ad_zwol(i)=0.0_r8
473 END DO
474 DO j=jmins,jmaxs
475 DO i=imins,imaxs
476 ad_hlv(i,j)=0.0_r8
477 ad_lheat(i,j)=0.0_r8
478 ad_lrad(i,j)=0.0_r8
479 ad_sheat(i,j)=0.0_r8
480 ad_taux(i,j)=0.0_r8
481 ad_tauy(i,j)=0.0_r8
482 END DO
483 END DO
484!
485 twopi_inv=0.5_r8/pi
486!
487!-----------------------------------------------------------------------
488! Exchange boundary data.
489!-----------------------------------------------------------------------
490!
491# ifdef DISTRIBUTE
492!^ CALL mp_exchange2d (ng, tile, iTLM, 2, &
493!^ & LBi, UBi, LBj, UBj, &
494!^ & NghostPoints, &
495!^ & EWperiodic(ng), NSperiodic(ng), &
496!^ & tl_sustr, tl_svstr)
497!^
498 CALL ad_mp_exchange2d (ng, tile, iadm, 2, &
499 & lbi, ubi, lbj, ubj, &
500 & nghostpoints, &
501 & ewperiodic(ng), nsperiodic(ng), &
502 & ad_sustr, ad_svstr)
503# ifdef EMINUSP
504!^ CALL mp_exchange2d (ng, tile, iTLM, 2, &
505!^ & LBi, UBi, LBj, UBj, &
506!^ & NghostPoints, &
507!^ & EWperiodic(ng), NSperiodic(ng), &
508!^ & tl_evap, &
509!^ & tl_stflx(:,:,isalt))
510!^
511 CALL ad_mp_exchange2d (ng, tile, iadm, 2, &
512 & lbi, ubi, lbj, ubj, &
513 & nghostpoints, &
514 & ewperiodic(ng), nsperiodic(ng), &
515 & ad_evap, &
516 & ad_stflx(:,:,isalt))
517# endif
518!^ CALL mp_exchange2d (ng, tile, iTLM, 4, &
519!^ & LBi, UBi, LBj, UBj, &
520!^ & NghostPoints, &
521!^ & EWperiodic(ng), NSperiodic(ng), &
522!^ & tl_lrflx, tl_lhflx, tl_shflx, &
523!^ & tl_stflx(:,:,itemp))
524!^
525 CALL ad_mp_exchange2d (ng, tile, iadm, 4, &
526 & lbi, ubi, lbj, ubj, &
527 & nghostpoints, &
528 & ewperiodic(ng), nsperiodic(ng), &
529 & ad_lrflx, ad_lhflx, ad_shflx, &
530 & ad_stflx(:,:,itemp))
531# endif
532 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
533!^ CALL exchange_v2d_tile (ng, tile, &
534!^ & LBi, UBi, LBj, UBj, &
535!^ & tl_svstr)
536!^
537 CALL ad_exchange_v2d_tile (ng, tile, &
538 & lbi, ubi, lbj, ubj, &
539 & ad_svstr)
540!^ CALL exchange_u2d_tile (ng, tile, &
541!^ & LBi, UBi, LBj, UBj, &
542!^ & tl_sustr)
543!^
544 CALL ad_exchange_u2d_tile (ng, tile, &
545 & lbi, ubi, lbj, ubj, &
546 & ad_sustr)
547# ifdef EMINUSP
548!^ CALL exchange_r2d_tile (ng, tile, &
549!^ & LBi, UBi, LBj, UBj, &
550!^ & tl_stflx(:,:,isalt))
551!^
552 CALL ad_exchange_r2d_tile (ng, tile, &
553 & lbi, ubi, lbj, ubj, &
554 & ad_stflx(:,:,isalt))
555!^ CALL exchange_r2d_tile (ng, tile, &
556!^ & LBi, UBi, LBj, UBj, &
557!^ & tl_evap)
558!^
559 CALL ad_exchange_r2d_tile (ng, tile, &
560 & lbi, ubi, lbj, ubj, &
561 & ad_evap)
562# endif
563!^ CALL exchange_r2d_tile (ng, tile, &
564!^ & LBi, UBi, LBj, UBj, &
565!^ & tl_stflx(:,:,itemp))
566!^
567 CALL ad_exchange_r2d_tile (ng, tile, &
568 & lbi, ubi, lbj, ubj, &
569 & ad_stflx(:,:,itemp))
570!^ CALL exchange_r2d_tile (ng, tile, &
571!^ & LBi, UBi, LBj, UBj, &
572!^ & tl_shflx)
573!^
574 CALL ad_exchange_r2d_tile (ng, tile, &
575 & lbi, ubi, lbj, ubj, &
576 & ad_shflx)
577!^ CALL exchange_r2d_tile (ng, tile, &
578!^ & LBi, UBi, LBj, UBj, &
579!^ & tl_lhflx)
580!^
581 CALL ad_exchange_r2d_tile (ng, tile, &
582 & lbi, ubi, lbj, ubj, &
583 & ad_lhflx)
584!^ CALL exchange_r2d_tile (ng, tile, &
585!^ & LBi, UBi, LBj, UBj, &
586!^ & tl_lrflx)
587!^
588 CALL ad_exchange_r2d_tile (ng, tile, &
589 & lbi, ubi, lbj, ubj, &
590 & ad_lrflx)
591 END IF
592!
593!=======================================================================
594! Compute adjoint surface net heat flux and surface wind stress.
595!=======================================================================
596!
597! Compute kinematic, surface, net heat flux (degC m/s). Notice that
598! the signs of latent and sensible fluxes are reversed because fluxes
599! calculated from the bulk formulations above are positive out of the
600! ocean.
601!
602! For EMINUSP option, EVAP = LHeat (W/m2) / Hlv (J/kg) = kg/m2/s
603! PREC = rain = kg/m2/s
604!
605! To convert these rates to m/s divide by freshwater density, rhow.
606!
607! Note that when the air is undersaturated in water vapor (Q < Qsea)
608! the model will evaporate and LHeat > 0:
609!
610! LHeat positive out of the ocean
611! evap positive out of the ocean
612!
613! Note that if evaporating, the salt flux is positive
614! and if raining, the salt flux is negative
615!
616! Note that fresh water flux is positive out of the ocean and the
617! salt flux (stflx(isalt)) is positive into the ocean. It is converted
618! to (psu m/s) for stflx(isalt) in "set_vbc.F".
619!
620
621
622!
623! Compute adjoint of kinematic, surface wind stress (m2/s2).
624!
625 cff=0.5_r8/rho0
626 DO j=jstrr,jendr
627 DO i=istr,iendr
628# ifdef MASKING
629!^ tl_sustr(i,j)=tl_sustr(i,j)*umask(i,j)
630 ad_sustr(i,j)=ad_sustr(i,j)*umask(i,j)
631# endif
632!^ tl_sustr(i,j)=cff*(tl_Taux(i-1,j)+tl_Taux(i,j))
633 adfac=cff*ad_sustr(i,j)
634 ad_taux(i-1,j)=ad_taux(i-1,j)+adfac
635 ad_taux(i ,j)=ad_taux(i ,j)+adfac
636# if !(defined AD_SENSITIVITY || defined ADJUST_WSTRESS || \
637 defined i4dvar_ana_sensitivity || defined opt_observations || \
638 defined so_semi)
639 ad_sustr(i,j)=0.0_r8
640# endif
641 END DO
642 END DO
643 DO j=jstr,jendr
644 DO i=istrr,iendr
645# ifdef MASKING
646!^ tl_svstr(i,j)=tl_svstr(i,j)*vmask(i,j)
647 ad_svstr(i,j)=ad_svstr(i,j)*vmask(i,j)
648# endif
649!^ tl_svstr(i,j)=cff*(tl_Tauy(i,j-1)+tl_Tauy(i,j))
650 adfac=cff*ad_svstr(i,j)
651 ad_tauy(i,j-1)=ad_tauy(i,j-1)+adfac
652 ad_tauy(i,j )=ad_tauy(i,j )+adfac
653# if !(defined AD_SENSITIVITY || defined ADJUST_WSTRESS || \
654 defined i4dvar_ana_sensitivity || defined opt_observations || \
655 defined so_semi)
656 ad_svstr(i,j)=0.0_r8
657# endif
658 END DO
659 END DO
660
661!
662! Compute latent heat of vaporization (J/kg) at sea surface, Hlv.
663!
664 DO j=jstr-1,jendr
665 DO i=istr-1,iendr
666 tseac(i)=t(i,j,n(ng),nrhs,itemp)
667 hlv(i,j)=(2.501_r8-0.00237_r8*tseac(i))*1.0e+6_r8
668 END DO
669 END DO
670
671 hscale=1.0_r8/(rho0*cp)
672 cff=1.0_r8/rhow
673 DO j=jstrr,jendr
674 DO i=istrr,iendr
675# ifdef MASKING
676# ifdef EMINUSP
677!^ tl_evap(i,j)=tl_evap(i,j)*rmask(i,j)
678 ad_evap(i,j)=ad_evap(i,j)*rmask(i,j)
679!^ tl_stflx(i,j,isalt)=tl_stflx(i,j,isalt)*rmask(i,j)
680 ad_stflx(i,j,isalt)=ad_stflx(i,j,isalt)*rmask(i,j)
681# endif
682!^ tl_stflx(i,j,itemp)=tl_stflx(i,j,itemp)*rmask(i,j)
683 ad_stflx(i,j,itemp)=ad_stflx(i,j,itemp)*rmask(i,j)
684# endif
685# ifdef EMINUSP
686!^ tl_stflx(i,j,isalt)=cff*tl_evap(i,j)
687 ad_evap(i,j)=ad_evap(i,j)+cff*ad_stflx(i,j,isalt)
688# if !(defined AD_SENSITIVITY || defined ADJUST_STFLUX || \
689 defined i4dvar_ana_sensitivity || defined opt_observations || \
690 defined so_semi)
691 ad_stflx(i,j,isalt)=0.0_r8
692# endif
693!^ tl_evap(i,j)=tl_LHeat(i,j)/Hlv(i,j)- &
694!^ & tl_Hlv(i,j)*evap(i,j)/Hlv(i,j)
695 adfac=ad_evap(i,j)/hlv(i,j)
696 ad_lheat(i,j)=ad_lheat(i,j)+adfac
697 ad_hlv(i,j)=ad_hlv(i,j)-adfac*evap(i,j)
698# if !(defined AD_SENSITIVITY || defind I4DVAR_ANA_SENSITIVITY || \
699 defined opt_observations || defined so_semi)
700 ad_evap(i,j)=0.0_r8
701# endif
702# endif
703!^ tl_stflx(i,j,itemp)=(tl_lrflx(i,j)+ &
704!^ & tl_lhflx(i,j)+tl_shflx(i,j))
705 ad_lrflx(i,j)=ad_lrflx(i,j)+ad_stflx(i,j,itemp)
706 ad_lhflx(i,j)=ad_lhflx(i,j)+ad_stflx(i,j,itemp)
707 ad_shflx(i,j)=ad_shflx(i,j)+ad_stflx(i,j,itemp)
708# if !(defined AD_SENSITIVITY || defined I4DVAR_ANA_SENSITIVITY || \
709 defined opt_observations || defined so_semi)
710 ad_stflx(i,j,itemp)=0.0_r8
711# endif
712!^ tl_shflx(i,j)=-tl_SHeat(i,j)*Hscale
713 ad_sheat(i,j)=ad_sheat(i,j)-ad_shflx(i,j)*hscale
714# if !(defined AD_SENSITIVITY || defined I4DVAR_ANA_SENSITIVITY || \
715 defined opt_observations || defined so_semi)
716 ad_shflx(i,j)=0.0_r8
717# endif
718!^ tl_lhflx(i,j)=-tl_LHeat(i,j)*Hscale
719 ad_lheat(i,j)=ad_lheat(i,j)-ad_lhflx(i,j)*hscale
720# if !(defined AD_SENSITIVITY || defined I4DVAR_ANA_SENSITIVITY || \
721 defined opt_observations || defined so_semi)
722 ad_lhflx(i,j)=0.0_r8
723# endif
724!^ tl_lrflx(i,j)=tl_LRad(i,j)*Hscale
725 ad_lrad(i,j)=ad_lrad(i,j)+ad_lrflx(i,j)*hscale
726 ad_lrflx(i,j)=0.0_r8
727
728 END DO
729 END DO
730!
731!=======================================================================
732! Adjoint Atmosphere-Ocean bulk fluxes parameterization.
733!=======================================================================
734!
735 hscale=rho0*cp
736 DO j=jstr-1,jendr
737!
738! Compute Basic State quantities. required
739!
740 DO i=istr-1,iendr
741!
742! Input bulk parameterization fields.
743!
744 wmag(i)=sqrt(uwind(i,j)*uwind(i,j)+vwind(i,j)*vwind(i,j))
745 pairm=pair(i,j)
746 tairc(i)=tair(i,j)
747 tairk(i)=tairc(i)+273.16_r8
748 tseac(i)=t(i,j,n(ng),nrhs,itemp)
749 tseak(i)=tseac(i)+273.16_r8
750 rhosea(i)=rho(i,j,n(ng))+1000.0_r8
751 rh=hair(i,j)
752 srad(i,j)=srflx(i,j)*hscale
753 tcff(i)=alpha(i,j)
754 scff(i)=beta(i,j)
755!
756! Initialize.
757!
758 deltc(i)=0.0_r8
759 delqc(i)=0.0_r8
760 lheat(i,j)=lhflx(i,j)*hscale
761 sheat(i,j)=shflx(i,j)*hscale
762 taur=0.0_r8
763 taux(i,j)=0.0_r8
764 tauy(i,j)=0.0_r8
765!
766!-----------------------------------------------------------------------
767! Compute net longwave radiation (W/m2), LRad.
768!-----------------------------------------------------------------------
769
770# if defined LONGWAVE
771!
772! Use Berliand (1952) formula to calculate net longwave radiation.
773! The equation for saturation vapor pressure is from Gill (Atmosphere-
774! Ocean Dynamics, pp 606). Here the coefficient in the cloud term
775! is assumed constant, but it is a function of latitude varying from
776! 1.0 at poles to 0.5 at the equator).
777!
778 cff=(0.7859_r8+0.03477_r8*tairc(i))/ &
779 & (1.0_r8+0.00412_r8*tairc(i))
780 e_sat=10.0_r8**cff ! saturation vapor pressure (hPa or mbar)
781 vap_p=e_sat*rh ! water vapor pressure (hPa or mbar)
782 cff2=tairk(i)*tairk(i)*tairk(i)
783 cff1=cff2*tairk(i)
784 lrad(i,j)=-emmiss*stefbo* &
785 & (cff1*(0.39_r8-0.05_r8*sqrt(vap_p))* &
786 & (1.0_r8-0.6823_r8*cloud(i,j)*cloud(i,j))+ &
787 & cff2*4.0_r8*(tseak(i)-tairk(i)))
788
789# elif defined LONGWAVE_OUT
790!
791! Treat input longwave data as downwelling radiation only and add
792! outgoing IR from model sea surface temperature.
793!
794 lrad(i,j)=lrflx(i,j)*hscale- &
795 & emmiss*stefbo*tseak(i)*tseak(i)*tseak(i)*tseak(i)
796
797# else
798 lrad(i,j)=lrflx(i,j)*hscale
799# endif
800# ifdef MASKING
801 lrad(i,j)=lrad(i,j)*rmask(i,j)
802# endif
803!
804!-----------------------------------------------------------------------
805! Compute specific humidities (kg/kg).
806!
807! note that Qair(i) is the saturation specific humidity at Tair
808! Q(i) is the actual specific humidity
809! Qsea(i) is the saturation specific humidity at Tsea
810!
811! Saturation vapor pressure in mb is first computed and then
812! converted to specific humidity in kg/kg
813!
814! The saturation vapor pressure is computed from Teten formula
815! using the approach of Buck (1981):
816!
817! Esat(mb) = (1.0007_r8+3.46E-6_r8*PairM(mb))*6.1121_r8*
818! EXP(17.502_r8*TairC(C)/(240.97_r8+TairC(C)))
819!
820! The ambient vapor is found from the definition of the
821! Relative humidity:
822!
823! RH = W/Ws*100 ~ E/Esat*100 E = RH/100*Esat if RH is in %
824! E = RH*Esat if RH fractional
825!
826! The specific humidity is then found using the relationship:
827!
828! Q = 0.622 E/(P + (0.622-1)e)
829!
830! Q(kg/kg) = 0.62197_r8*(E(mb)/(PairM(mb)-0.378_r8*E(mb)))
831!
832!-----------------------------------------------------------------------
833!
834! Compute air saturation vapor pressure (mb), using Teten formula.
835!
836 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
837 & exp(17.502_r8*tairc(i)/(240.97_r8+tairc(i)))
838!
839! Compute specific humidity at Saturation, Qair (kg/kg).
840!
841 qair(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff))
842!
843! Compute specific humidity, Q (kg/kg).
844!
845 IF (rh.lt.2.0_r8) THEN !RH fraction
846 cff=cff*rh !Vapor pres (mb)
847 q(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff)) !Spec hum (kg/kg)
848 ELSE !RH input was actually specific humidity in g/kg
849 q(i)=rh/1000.0_r8 !Spec Hum (kg/kg)
850 END IF
851!
852! Compute water saturation vapor pressure (mb), using Teten formula.
853!
854 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
855 & exp(17.502_r8*tseac(i)/(240.97_r8+tseac(i)))
856!
857! Vapor Pressure reduced for salinity (Kraus and Businger, 1994, pp42).
858!
859 cff=cff*0.98_r8
860!
861! Compute Qsea (kg/kg) from vapor pressure.
862!
863 qsea(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff))
864!
865!-----------------------------------------------------------------------
866! Compute Monin-Obukhov similarity parameters for wind (Wstar),
867! heat (Tstar), and moisture (Qstar), Liu et al. (1979).
868!-----------------------------------------------------------------------
869!
870! Moist air density (kg/m3).
871!
872 rhoair(i)=pairm*100.0_r8/(blk_rgas*tairk(i)* &
873 & (1.0_r8+0.61_r8*q(i)))
874!
875! Kinematic viscosity of dry air (m2/s), Andreas (1989).
876!
877 visair(i)=1.326e-5_r8* &
878 & (1.0_r8+tairc(i)*(6.542e-3_r8+tairc(i)* &
879 & (8.301e-6_r8-4.84e-9_r8*tairc(i))))
880!
881! Compute latent heat of vaporization (J/kg) at sea surface, Hlv.
882!
883 hlv(i,j)=(2.501_r8-0.00237_r8*tseac(i))*1.0e+6_r8
884!
885! Assume that wind is measured relative to sea surface and include
886! gustiness.
887!
888 wgus(i)=0.5_r8
889 delw(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
890 delq(i)=qsea(i)-q(i)
891 delt(i)=tseac(i)-tairc(i)
892!
893! Neutral coefficients.
894!
895 zow(i)=0.0001_r8
896 u10(i)=delw(i)*log(10.0_r8/zow(i))/log(blk_zw(ng)/zow(i))
897 wstar(i)=0.035_r8*u10(i)
898 zo10(i)=0.011_r8*wstar(i)*wstar(i)/g+ &
899 & 0.11_r8*visair(i)/wstar(i)
900 cd10(i)=(vonkar/log(10.0_r8/zo10(i)))**2
901 ch10(i)=0.00115_r8
902 ct10(i)=ch10(i)/sqrt(cd10(i))
903 zot10(i)=10.0_r8/exp(vonkar/ct10(i))
904 cd=(vonkar/log(blk_zw(ng)/zo10(i)))**2
905!
906! Compute Richardson number.
907!
908 ct(i)=vonkar/log(blk_zt(ng)/zot10(i)) ! T transfer coefficient
909 cc(i)=vonkar*ct(i)/cd
910 deltc(i)=0.0_r8
911# ifdef COOL_SKIN
912 deltc(i)=blk_dter
913# endif
914 ribcu(i)=-blk_zw(ng)/(blk_zabl*0.004_r8*blk_beta**3)
915 ri(i)=-g*blk_zw(ng)*((delt(i)-deltc(i))+ &
916 & 0.61_r8*tairk(i)*delq(i))/ &
917 & (tairk(i)*delw(i)*delw(i))
918 IF (ri(i).lt.0.0_r8) THEN
919 zetu(i)=cc(i)*ri(i)/(1.0_r8+ri(i)/ribcu(i)) ! Unstable
920 ELSE
921 zetu(i)=cc(i)*ri(i)/(1.0_r8+3.0_r8*ri(i)/cc(i)) ! Stable
922 END IF
923 l10(i)=blk_zw(ng)/zetu(i)
924!
925! First guesses for Monon-Obukhov similarity scales.
926!
927 wstar(i)=delw(i)*vonkar/(log(blk_zw(ng)/zo10(i))- &
928 & bulk_psiu(blk_zw(ng)/l10(i),pi))
929 tstar(i)=-(delt(i)-deltc(i))*vonkar/ &
930 & (log(blk_zt(ng)/zot10(i))- &
931 & bulk_psit(blk_zt(ng)/l10(i),pi))
932 qstar(i)=-(delq(i)-delqc(i))*vonkar/ &
933 & (log(blk_zq(ng)/zot10(i))- &
934 & bulk_psit(blk_zq(ng)/l10(i),pi))
935!
936! Modify Charnock for high wind speeds. The 0.125 factor below is for
937! 1.0/(18.0-10.0).
938!
939 IF (delw(i).gt.18.0_r8) THEN
940 charn(i)=0.018_r8
941 ELSE IF ((10.0_r8.lt.delw(i)).and.(delw(i).le.18.0_r8)) THEN
942 charn(i)=0.011_r8+ &
943 & 0.125_r8*(0.018_r8-0.011_r8)*(delw(i)-10._r8)
944 ELSE
945 charn(i)=0.011_r8
946 END IF
947# ifdef BBL_MODEL_NOT_YET
948 cwave(i)=g*pwave(i,j)*twopi_inv
949 lwave(i)=cwave(i)*pwave(i,j)
950# endif
951 END DO
952!
953! Iterate until convergence. It usually converges within 3 iterations.
954!
955 DO iter=1,itermax
956 DO i=istr-1,iendr
957# ifdef BBL_MODEL_NOT_YET
958!
959! Use wave info if we have it, two different options.
960!
961# ifdef WIND_WAVES
962 zow(i)=(25._r8/pi)*lwave(i)*(wstar(i)/cwave(i))**4.5+ &
963 & 0.11_r8*visair(i)/(wstar(i)+eps)
964# else
965 zow(i)=1200._r8*awave(i,j)*(awave(i,j)/lwave(i))**4.5+ &
966 & 0.11_r8*visair(i)/(wstar(i)+eps)
967# endif
968# else
969 zow(i)=charn(i)*wstar(i)*wstar(i)/g+ &
970 & 0.11_r8*visair(i)/(wstar(i)+eps)
971# endif
972 rr(i)=zow(i)*wstar(i)/visair(i)
973!
974! Compute Monin-Obukhov stability parameter, Z/L.
975!
976 zoq(i)=min(1.15e-4_r8,5.5e-5_r8/rr(i)**0.6)
977 zot(i)=zoq(i)
978 zol(i)=vonkar*g*blk_zw(ng)* &
979 & (tstar(i)*(1.0_r8+0.61_r8*q(i))+ &
980 & 0.61_r8*tairk(i)*qstar(i))/ &
981 & (tairk(i)*wstar(i)*wstar(i)* &
982 & (1.0_r8+0.61_r8*q(i))+eps)
983 l(i)=blk_zw(ng)/(zol(i)+eps)
984!
985! Evaluate stability functions at Z/L.
986!
987 wpsi(i)=bulk_psiu(zol(i),pi)
988 tpsi(i)=bulk_psit(blk_zt(ng)/l(i),pi)
989 qpsi(i)=bulk_psit(blk_zq(ng)/l(i),pi)
990# ifdef COOL_SKIN
991 cwet(i)=0.622_r8*hlv(i,j)*qsea(i)/ &
992 & (blk_rgas*tseak(i)*tseak(i))
993 delqc(i)=cwet(i)*deltc(i)
994# endif
995!
996! Compute wind scaling parameters, Wstar.
997!
998 wstar(i)=max(eps,delw(i)*vonkar/ &
999 & (log(blk_zw(ng)/zow(i))-wpsi(i)))
1000 tstar(i)=-(delt(i)-deltc(i))*vonkar/ &
1001 & (log(blk_zt(ng)/zot(i))-tpsi(i))
1002 qstar(i)=-(delq(i)-delqc(i))*vonkar/ &
1003 & (log(blk_zq(ng)/zoq(i))-qpsi(i))
1004!
1005! Compute gustiness in wind speed.
1006!
1007 bf=-g/tairk(i)* &
1008 & wstar(i)*(tstar(i)+0.61_r8*tairk(i)*qstar(i))
1009 IF (bf.gt.0.0_r8) THEN
1010 wgus(i)=blk_beta*(bf*blk_zabl)**r3
1011 ELSE
1012 wgus(i)=0.2_r8
1013 END IF
1014 delw(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
1015# ifdef COOL_SKIN
1016!
1017!-----------------------------------------------------------------------
1018! Cool Skin correction.
1019!-----------------------------------------------------------------------
1020!
1021! Cool skin correction constants. Clam: part of Saunders constant
1022! lambda; Cwet: slope of saturation vapor.
1023!
1024 clam=16.0_r8*g*blk_cpw*(rhosea(i)*blk_visw)**3/ &
1025 & (blk_tcw*blk_tcw*rhoair(i)*rhoair(i))
1026!
1027! Set initial guesses for cool-skin layer thickness (Hcool).
1028!
1029 hcool=0.001_r8
1030!
1031! Backgound sensible and latent heat.
1032!
1033 hsb=-rhoair(i)*blk_cpa*wstar(i)*tstar(i)
1034 hlb=-rhoair(i)*hlv(i,j)*wstar(i)*qstar(i)
1035!
1036! Mean absoption in cool-skin layer.
1037!
1038 fc=0.065_r8+11.0_r8*hcool- &
1039 & (1.0_r8-exp(-hcool*1250.0_r8))*6.6e-5_r8/hcool
1040!
1041! Total cooling at the interface.
1042!
1043 qcool=lrad(i,j)+hsb+hlb-srad(i,j)*fc
1044 qbouy=tcff(i)*qcool+scff(i)*hlb*blk_cpw/hlv(i,j)
1045!
1046! Compute temperature and moisture change.
1047!
1048 IF ((qcool.gt.0.0_r8).and.(qbouy.gt.0.0_r8)) THEN
1049 lambd=6.0_r8/ &
1050 & (1.0_r8+(clam*qbouy/(wstar(i)+eps)**4)**0.75_r8)**r3
1051 hcool=lambd*blk_visw/(sqrt(rhoair(i)/rhosea(i))* &
1052 & wstar(i)+eps)
1053 deltc(i)=qcool*hcool/blk_tcw
1054 ELSE
1055 deltc(i)=0.0_r8
1056 END IF
1057 delqc(i)=cwet(i)*deltc(i)
1058# endif
1059 END DO
1060 END DO
1061!
1062
1063!
1064!-----------------------------------------------------------------------
1065! Compute Adjoint of Atmosphere/Ocean fluxes.
1066!-----------------------------------------------------------------------
1067!
1068 DO i=istr-1,iendr
1069!
1070! Compute transfer coefficients for momentum (Cd).
1071!
1072 wspeed=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
1073 cd=wstar(i)*wstar(i)/(wspeed*wspeed+eps)
1074!
1075! Compute turbulent sensible heat flux (W/m2), Hs.
1076!
1077 hs=-blk_cpa*rhoair(i)*wstar(i)*tstar(i)
1078!
1079! Compute sensible heat flux (W/m2) due to rainfall (kg/m2/s), Hsr.
1080!
1081 diffw=2.11e-5_r8*(tairk(i)/273.16_r8)**1.94_r8
1082 diffh=0.02411_r8*(1.0_r8+tairc(i)* &
1083 & (3.309e-3_r8-1.44e-6_r8*tairc(i)))/ &
1084 & (rhoair(i)*blk_cpa)
1085 cff=qair(i)*hlv(i,j)/(blk_rgas*tairk(i)*tairk(i))
1086 wet_bulb=1.0_r8/(1.0_r8+0.622_r8*(cff*hlv(i,j)*diffw)/ &
1087 & (blk_cpa*diffh))
1088 hsr=rain(i,j)*wet_bulb*blk_cpw* &
1089 & ((tseac(i)-tairc(i))+(qsea(i)-q(i))*hlv(i,j)/blk_cpa)
1090 sheat(i,j)=(hs+hsr)
1091# ifdef MASKING
1092 sheat(i,j)=sheat(i,j)*rmask(i,j)
1093# endif
1094!
1095! Compute turbulent latent heat flux (W/m2), Hl.
1096!
1097 hl=-hlv(i,j)*rhoair(i)*wstar(i)*qstar(i)
1098!
1099! Compute Webb correction (Webb effect) to latent heat flux, Hlw.
1100!
1101 upvel=-1.61_r8*wstar(i)*qstar(i)- &
1102 & (1.0_r8+1.61_r8*q(i))*wstar(i)*tstar(i)/tairk(i)
1103 hlw=rhoair(i)*hlv(i,j)*upvel*q(i)
1104 lheat(i,j)=(hl+hlw)
1105# ifdef MASKING
1106 lheat(i,j)=lheat(i,j)*rmask(i,j)
1107# endif
1108!
1109! Adjoint wind stress components (N/m2), Tau.
1110!
1111
1112# ifdef MASKING
1113!^ tl_Tauy(i,j)=tl_Tauy(i,j)*rmask(i,j)
1114 ad_tauy(i,j)=ad_tauy(i,j)*rmask(i,j)
1115# endif
1116!^ tl_Tauy(i,j)=tl_cff*Vwind(i,j)
1117 ad_cff=ad_cff+ad_tauy(i,j)*vwind(i,j)
1118 ad_tauy(i,j)=0.0_r8
1119# ifdef MASKING
1120!^ tl_Taux(i,j)=tl_Taux(i,j)*rmask(i,j)
1121 ad_taux(i,j)=ad_taux(i,j)*rmask(i,j)
1122# endif
1123!^ tl_Taux(i,j)=tl_cff*Uwind(i,j)
1124 ad_cff=ad_cff+ad_taux(i,j)*uwind(i,j)
1125 ad_taux(i,j)=0.0_r8
1126
1127!^ tl_cff=rhoAir(i)*(tl_Cd*Wspeed+Cd*tl_Wspeed)
1128 adfac=rhoair(i)*ad_cff
1129 ad_cd=ad_cd+wspeed*adfac
1130 ad_wspeed=ad_wspeed+cd*adfac
1131 ad_cff=0.0_r8
1132!
1133! Adjoint turbulent latent heat flux (W/m2), Hl.
1134!
1135
1136
1137# ifdef MASKING
1138!^ tl_LHeat(i,j)=tl_LHeat(i,j)*rmask(i,j)
1139 ad_lheat(i,j)=ad_lheat(i,j)*rmask(i,j)
1140# endif
1141!^ tl_LHeat(i,j)=(tl_Hl+tl_Hlw)
1142 ad_hl=ad_hl+ad_lheat(i,j)
1143 ad_hlw=ad_hlw+ad_lheat(i,j)
1144 ad_lheat(i,j)=0.0_r8
1145
1146!^ tl_Hlw=rhoAir(i)*Q(i)*(tl_Hlv(i,j)*upvel+ &
1147!^ & Hlv(i,j)*tl_upvel)
1148 adfac=rhoair(i)*q(i)*ad_hlw
1149 ad_hlv(i,j)=ad_hlv(i,j)+upvel*adfac
1150 ad_upvel=ad_upvel+hlv(i,j)*adfac
1151 ad_hlw=0.0_r8
1152
1153!^ tl_upvel=-1.61_r8* &
1154!^ & (tl_Wstar(i)*Qstar(i)+Wstar(i)*tl_Qstar(i))- &
1155!^ & (1.0_r8+1.61_r8*Q(i))*(tl_Wstar(i)*Tstar(i)+ &
1156!^ & Wstar(i)*tl_Tstar(i))/TairK(i)
1157 adfac=-1.61_r8*ad_upvel
1158 adfac1=-(1.0_r8+1.61_r8*q(i))*ad_upvel/tairk(i)
1159 ad_wstar(i)=ad_wstar(i)+qstar(i)*adfac+tstar(i)*adfac1
1160 ad_qstar(i)=ad_qstar(i)+wstar(i)*adfac
1161 ad_tstar(i)=ad_tstar(i)+wstar(i)*adfac1
1162 ad_upvel=0.0_r8
1163
1164!^ tl_Hl=-tl_Hlv(i,j)*rhoAir(i)*Wstar(i)*Qstar(i)- &
1165!^ & Hlv(i,j)*rhoAir(i)*(tl_Wstar(i)*Qstar(i)+ &
1166!^ & Wstar(i)*tl_Qstar(i) )
1167 adfac=hlv(i,j)*rhoair(i)*ad_hl
1168 ad_hlv(i,j)=ad_hlv(i,j)-rhoair(i)*wstar(i)*qstar(i)*ad_hl
1169 ad_wstar(i)=ad_wstar(i)-qstar(i)*adfac
1170 ad_qstar(i)=ad_qstar(i)-wstar(i)*adfac
1171 ad_hl=0.0_r8
1172!
1173! Adjoint sensible heat flux (W/m2) due to rainfall (kg/m2/s), Hsr.
1174!
1175
1176
1177# ifdef MASKING
1178!^ tl_SHeat(i,j)=tl_SHeat(i,j)*rmask(i,j)
1179 ad_sheat(i,j)=ad_sheat(i,j)*rmask(i,j)
1180# endif
1181
1182!^ tl_SHeat(i,j)=(tl_Hs+tl_Hsr)
1183 ad_hs=ad_hs+ad_sheat(i,j)
1184 ad_hsr=ad_hsr+ad_sheat(i,j)
1185 ad_sheat(i,j)=0.0_r8
1186
1187!^ tl_Hsr=Hsr*tl_wet_bulb/wet_bulb+ &
1188!^ & rain(i,j)*wet_bulb*blk_Cpw* &
1189!^ & (tl_TseaC(i)+tl_Qsea(i)*Hlv(i,j)/blk_Cpa+ &
1190!^ & (Qsea(i)-Q(i))*tl_Hlv(i,j)/blk_Cpa)
1191 adfac=rain(i,j)*wet_bulb*blk_cpw*ad_hsr
1192 adfac1=adfac/blk_cpa
1193 ad_wet_bulb=ad_wet_bulb+hsr*ad_hsr/wet_bulb
1194 ad_tseac(i)=ad_tseac(i)+adfac
1195 ad_qsea(i)=ad_qsea(i)+hlv(i,j)*adfac1
1196 ad_hlv(i,j)=ad_hlv(i,j)+(qsea(i)-q(i))*adfac1
1197 ad_hsr=0.0_r8
1198
1199!^ tl_wet_bulb=-tl_fac*wet_bulb*wet_bulb
1200 ad_fac=-wet_bulb*wet_bulb*ad_wet_bulb
1201 ad_wet_bulb=0.0_r8
1202
1203!^ tl_fac=0.622_r8*diffw*(tl_cff*Hlv(i,j)+cff*tl_Hlv(i,j))/ &
1204!^ & (blk_Cpa*diffh)
1205 adfac=0.622_r8*diffw*ad_fac/(blk_cpa*diffh)
1206 ad_cff=ad_cff+hlv(i,j)*adfac
1207 ad_hlv(i,j)=ad_hlv(i,j)+cff*adfac
1208 ad_fac=0.0_r8
1209
1210!^ tl_cff=Qair(i)*tl_Hlv(i,j)/(blk_Rgas*TairK(i)*TairK(i))
1211 ad_hlv(i,j)=ad_hlv(i,j)+qair(i)*ad_cff/ &
1212 & (blk_rgas*tairk(i)*tairk(i))
1213 ad_cff=0.0_r8
1214!
1215! Adjoint turbulent sensible heat flux (W/m2), Hs.
1216!
1217
1218!^ tl_Hs=-blk_Cpa*rhoAir(i)*(tl_Wstar(i)*Tstar(i)+ &
1219!^ & Wstar(i)*tl_Tstar(i))
1220 adfac=-blk_cpa*rhoair(i)*ad_hs
1221 ad_wstar(i)=ad_wstar(i)+tstar(i)*adfac
1222 ad_tstar(i)=ad_tstar(i)+wstar(i)*adfac
1223 ad_hs=0.0_r8
1224
1225!
1226! Adjoint transfer coefficients for momentum (Cd).
1227!
1228!^ tl_Cd=2.0_r8*(tl_Wstar(i)*Wstar(i)/(Wspeed*Wspeed+eps) &
1229!^ & -tl_Wspeed*Wspeed*Cd/(Wspeed*Wspeed+eps))
1230 adfac=2.0_r8*ad_cd/(wspeed*wspeed+eps)
1231 ad_wstar(i)=ad_wstar(i)+wstar(i)*adfac
1232 ad_wspeed=ad_wspeed-wspeed*cd*adfac
1233 ad_cd=0.0_r8
1234
1235 IF (wspeed.ne.0.0_r8) THEN
1236!^ tl_Wspeed=tl_Wgus(i)*Wgus(i)/Wspeed
1237 ad_wgus(i)=ad_wgus(i)+wgus(i)*ad_wspeed/wspeed
1238 ad_wspeed=0.0_r8
1239 ELSE
1240!^ tl_Wspeed=0.0_r8
1241 ad_wspeed=0.0_r8
1242 END IF
1243
1244 END DO
1245!
1246! Adjoint Iteration until convergence. It usually converges within
1247! 3 iterations.
1248!
1249 DO iter=itermax,1,-1
1250!
1251! Recompute appropriate basic state variables again prior to
1252! each adjoint iteration.
1253!
1254 DO i=istr-1,iendr
1255!
1256! Input bulk parameterization fields.
1257!
1258 wmag(i)=sqrt(uwind(i,j)*uwind(i,j)+vwind(i,j)*vwind(i,j))
1259 pairm=pair(i,j)
1260 tairc(i)=tair(i,j)
1261 tairk(i)=tairc(i)+273.16_r8
1262 tseac(i)=t(i,j,n(ng),nrhs,itemp)
1263 tseak(i)=tseac(i)+273.16_r8
1264 rhosea(i)=rho(i,j,n(ng))+1000.0_r8
1265 rh=hair(i,j)
1266 srad(i,j)=srflx(i,j)*hscale
1267 tcff(i)=alpha(i,j)
1268 scff(i)=beta(i,j)
1269!
1270! Initialize.
1271!
1272 deltc(i)=0.0_r8
1273 delqc(i)=0.0_r8
1274 lheat(i,j)=lhflx(i,j)*hscale
1275 sheat(i,j)=shflx(i,j)*hscale
1276 taur=0.0_r8
1277 taux(i,j)=0.0_r8
1278 tauy(i,j)=0.0_r8
1279!
1280!-----------------------------------------------------------------------
1281! Compute net longwave radiation (W/m2), LRad.
1282!-----------------------------------------------------------------------
1283
1284# if defined LONGWAVE
1285!
1286! Use Berliand (1952) formula to calculate net longwave radiation.
1287! The equation for saturation vapor pressure is from Gill (Atmosphere-
1288! Ocean Dynamics, pp 606). Here the coefficient in the cloud term
1289! is assumed constant, but it is a function of latitude varying from
1290! 1.0 at poles to 0.5 at the equator).
1291!
1292 cff=(0.7859_r8+0.03477_r8*tairc(i))/ &
1293 & (1.0_r8+0.00412_r8*tairc(i))
1294 e_sat=10.0_r8**cff ! saturation vapor pressure (hPa or mbar)
1295 vap_p=e_sat*rh ! water vapor pressure (hPa or mbar)
1296 cff2=tairk(i)*tairk(i)*tairk(i)
1297 cff1=cff2*tairk(i)
1298 lrad(i,j)=-emmiss*stefbo* &
1299 & (cff1*(0.39_r8-0.05_r8*sqrt(vap_p))* &
1300 & (1.0_r8-0.6823_r8*cloud(i,j)*cloud(i,j))+ &
1301 & cff2*4.0_r8*(tseak(i)-tairk(i)))
1302
1303# elif defined LONGWAVE_OUT
1304!
1305! Treat input longwave data as downwelling radiation only and add
1306! outgoing IR from model sea surface temperature.
1307!
1308 lrad(i,j)=lrflx(i,j)*hscale- &
1309 & emmiss*stefbo*tseak(i)*tseak(i)*tseak(i)*tseak(i)
1310
1311# else
1312 lrad(i,j)=lrflx(i,j)*hscale
1313# endif
1314# ifdef MASKING
1315 lrad(i,j)=lrad(i,j)*rmask(i,j)
1316# endif
1317!
1318!-----------------------------------------------------------------------
1319! Compute specific humidities (kg/kg).
1320!
1321! note that Qair(i) is the saturation specific humidity at Tair
1322! Q(i) is the actual specific humidity
1323! Qsea(i) is the saturation specific humidity at Tsea
1324!
1325! Saturation vapor pressure in mb is first computed and then
1326! converted to specific humidity in kg/kg
1327!
1328! The saturation vapor pressure is computed from Teten formula
1329! using the approach of Buck (1981):
1330!
1331! Esat(mb) = (1.0007_r8+3.46E-6_r8*PairM(mb))*6.1121_r8*
1332! EXP(17.502_r8*TairC(C)/(240.97_r8+TairC(C)))
1333!
1334! The ambient vapor is found from the definition of the
1335! Relative humidity:
1336!
1337! RH = W/Ws*100 ~ E/Esat*100 E = RH/100*Esat if RH is in %
1338! E = RH*Esat if RH fractional
1339!
1340! The specific humidity is then found using the relationship:
1341!
1342! Q = 0.622 E/(P + (0.622-1)e)
1343!
1344! Q(kg/kg) = 0.62197_r8*(E(mb)/(PairM(mb)-0.378_r8*E(mb)))
1345!
1346!-----------------------------------------------------------------------
1347!
1348! Compute air saturation vapor pressure (mb), using Teten formula.
1349!
1350 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
1351 & exp(17.502_r8*tairc(i)/(240.97_r8+tairc(i)))
1352!
1353! Compute specific humidity at Saturation, Qair (kg/kg).
1354!
1355 qair(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff))
1356!
1357! Compute specific humidity, Q (kg/kg).
1358!
1359 IF (rh.lt.2.0_r8) THEN !RH fraction
1360 cff=cff*rh !Vapor pres (mb)
1361 q(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff)) !Spec hum (kg/kg)
1362 ELSE !RH input was actually specific humidity in g/kg
1363 q(i)=rh/1000.0_r8 !Spec Hum (kg/kg)
1364 END IF
1365!
1366! Compute water saturation vapor pressure (mb), using Teten formula.
1367!
1368 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
1369 & exp(17.502_r8*tseac(i)/(240.97_r8+tseac(i)))
1370!
1371! Vapor Pressure reduced for salinity (Kraus and Businger, 1994, pp42).
1372!
1373 cff=cff*0.98_r8
1374!
1375! Compute Qsea (kg/kg) from vapor pressure.
1376!
1377 qsea(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff))
1378!
1379!-----------------------------------------------------------------------
1380! Compute Monin-Obukhov similarity parameters for wind (Wstar),
1381! heat (Tstar), and moisture (Qstar), Liu et al. (1979).
1382!-----------------------------------------------------------------------
1383!
1384! Moist air density (kg/m3).
1385!
1386 rhoair(i)=pairm*100.0_r8/(blk_rgas*tairk(i)* &
1387 & (1.0_r8+0.61_r8*q(i)))
1388!
1389! Kinematic viscosity of dry air (m2/s), Andreas (1989).
1390!
1391 visair(i)=1.326e-5_r8* &
1392 & (1.0_r8+tairc(i)*(6.542e-3_r8+tairc(i)* &
1393 & (8.301e-6_r8-4.84e-9_r8*tairc(i))))
1394!
1395! Compute latent heat of vaporization (J/kg) at sea surface, Hlv.
1396!
1397 hlv(i,j)=(2.501_r8-0.00237_r8*tseac(i))*1.0e+6_r8
1398!
1399! Assume that wind is measured relative to sea surface and include
1400! gustiness.
1401!
1402 wgus(i)=0.5_r8
1403 delw(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
1404 delq(i)=qsea(i)-q(i)
1405 delt(i)=tseac(i)-tairc(i)
1406!
1407! Neutral coefficients.
1408!
1409 zow(i)=0.0001_r8
1410 u10(i)=delw(i)*log(10.0_r8/zow(i))/log(blk_zw(ng)/zow(i))
1411 wstar(i)=0.035_r8*u10(i)
1412 zo10(i)=0.011_r8*wstar(i)*wstar(i)/g+ &
1413 & 0.11_r8*visair(i)/wstar(i)
1414 cd10(i)=(vonkar/log(10.0_r8/zo10(i)))**2
1415 ch10(i)=0.00115_r8
1416 ct10(i)=ch10(i)/sqrt(cd10(i))
1417 zot10(i)=10.0_r8/exp(vonkar/ct10(i))
1418 cd=(vonkar/log(blk_zw(ng)/zo10(i)))**2
1419!
1420! Compute Richardson number.
1421!
1422 ct(i)=vonkar/log(blk_zt(ng)/zot10(i)) ! T transfer coefficient
1423 cc(i)=vonkar*ct(i)/cd
1424 deltc(i)=0.0_r8
1425# ifdef COOL_SKIN
1426 deltc(i)=blk_dter
1427# endif
1428 ribcu(i)=-blk_zw(ng)/(blk_zabl*0.004_r8*blk_beta**3)
1429 ri(i)=-g*blk_zw(ng)*((delt(i)-deltc(i))+ &
1430 & 0.61_r8*tairk(i)*delq(i))/ &
1431 & (tairk(i)*delw(i)*delw(i))
1432 IF (ri(i).lt.0.0_r8) THEN
1433 zetu(i)=cc(i)*ri(i)/(1.0_r8+ri(i)/ribcu(i)) ! Unstable
1434 ELSE
1435 zetu(i)=cc(i)*ri(i)/(1.0_r8+3.0_r8*ri(i)/cc(i)) ! Stable
1436 END IF
1437 l10(i)=blk_zw(ng)/zetu(i)
1438!
1439! First guesses for Monon-Obukhov similarity scales.
1440!
1441 wstar(i)=delw(i)*vonkar/(log(blk_zw(ng)/zo10(i))- &
1442 & bulk_psiu(blk_zw(ng)/l10(i),pi))
1443 tstar(i)=-(delt(i)-deltc(i))*vonkar/ &
1444 & (log(blk_zt(ng)/zot10(i))- &
1445 & bulk_psit(blk_zt(ng)/l10(i),pi))
1446 qstar(i)=-(delq(i)-delqc(i))*vonkar/ &
1447 & (log(blk_zq(ng)/zot10(i))- &
1448 & bulk_psit(blk_zq(ng)/l10(i),pi))
1449!
1450! Modify Charnock for high wind speeds. The 0.125 factor below is for
1451! 1.0/(18.0-10.0).
1452!
1453 IF (delw(i).gt.18.0_r8) THEN
1454 charn(i)=0.018_r8
1455 ELSE IF ((10.0_r8.lt.delw(i)).and.(delw(i).le.18.0_r8)) THEN
1456 charn(i)=0.011_r8+ &
1457 & 0.125_r8*(0.018_r8-0.011_r8)*(delw(i)-10._r8)
1458 ELSE
1459 charn(i)=0.011_r8
1460 END IF
1461# ifdef BBL_MODEL_NOT_YET
1462 cwave(i)=g*pwave(i,j)*twopi_inv
1463 lwave(i)=cwave(i)*pwave(i,j)
1464# endif
1465 END DO
1466!
1467! Iterate here to obtain the appropriate basic state variables.
1468!
1469 DO itera=1,iter-1
1470 DO i=istr-1,iendr
1471# ifdef BBL_MODEL_NOT_YET
1472!
1473! Use wave info if we have it, two different options.
1474!
1475# ifdef WIND_WAVES
1476 zow(i)=(25._r8/pi)*lwave(i)*(wstar(i)/cwave(i))**4.5+ &
1477 & 0.11_r8*visair(i)/(wstar(i)+eps)
1478# else
1479 zow(i)=1200._r8*awave(i,j)*(awave(i,j)/lwave(i))**4.5+ &
1480 & 0.11_r8*visair(i)/(wstar(i)+eps)
1481# endif
1482# else
1483 zow(i)=charn(i)*wstar(i)*wstar(i)/g+ &
1484 & 0.11_r8*visair(i)/(wstar(i)+eps)
1485# endif
1486 rr(i)=zow(i)*wstar(i)/visair(i)
1487!
1488! Compute Monin-Obukhov stability parameter, Z/L.
1489!
1490 zoq(i)=min(1.15e-4_r8,5.5e-5_r8/rr(i)**0.6)
1491 zot(i)=zoq(i)
1492 zol(i)=vonkar*g*blk_zw(ng)* &
1493 & (tstar(i)*(1.0_r8+0.61_r8*q(i))+ &
1494 & 0.61_r8*tairk(i)*qstar(i))/ &
1495 & (tairk(i)*wstar(i)*wstar(i)* &
1496 & (1.0_r8+0.61_r8*q(i))+eps)
1497 l(i)=blk_zw(ng)/(zol(i)+eps)
1498!
1499! Evaluate stability functions at Z/L.
1500!
1501 wpsi(i)=bulk_psiu(zol(i),pi)
1502 tpsi(i)=bulk_psit(blk_zt(ng)/l(i),pi)
1503 qpsi(i)=bulk_psit(blk_zq(ng)/l(i),pi)
1504# ifdef COOL_SKIN
1505 cwet(i)=0.622_r8*hlv(i,j)*qsea(i)/ &
1506 & (blk_rgas*tseak(i)*tseak(i))
1507 delqc(i)=cwet(i)*deltc(i)
1508# endif
1509!
1510! Compute wind scaling parameters, Wstar.
1511!
1512 wstar(i)=max(eps,delw(i)*vonkar/ &
1513 & (log(blk_zw(ng)/zow(i))-wpsi(i)))
1514 tstar(i)=-(delt(i)-deltc(i))*vonkar/ &
1515 & (log(blk_zt(ng)/zot(i))-tpsi(i))
1516 qstar(i)=-(delq(i)-delqc(i))*vonkar/ &
1517 & (log(blk_zq(ng)/zoq(i))-qpsi(i))
1518!
1519! Compute gustiness in wind speed.
1520!
1521 bf=-g/tairk(i)* &
1522 & wstar(i)*(tstar(i)+0.61_r8*tairk(i)*qstar(i))
1523 IF (bf.gt.0.0_r8) THEN
1524 wgus(i)=blk_beta*(bf*blk_zabl)**r3
1525 ELSE
1526 wgus(i)=0.2_r8
1527 END IF
1528 delw(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
1529# ifdef COOL_SKIN
1530!
1531!-----------------------------------------------------------------------
1532! Cool Skin correction.
1533!-----------------------------------------------------------------------
1534!
1535! Cool skin correction constants. Clam: part of Saunders constant
1536! lambda; Cwet: slope of saturation vapor.
1537!
1538 clam=16.0_r8*g*blk_cpw*(rhosea(i)*blk_visw)**3/ &
1539 & (blk_tcw*blk_tcw*rhoair(i)*rhoair(i))
1540!
1541! Set initial guesses for cool-skin layer thickness (Hcool).
1542!
1543 hcool=0.001_r8
1544!
1545! Backgound sensible and latent heat.
1546!
1547 hsb=-rhoair(i)*blk_cpa*wstar(i)*tstar(i)
1548 hlb=-rhoair(i)*hlv(i,j)*wstar(i)*qstar(i)
1549!
1550! Mean absoption in cool-skin layer.
1551!
1552 fc=0.065_r8+11.0_r8*hcool- &
1553 & (1.0_r8-exp(-hcool*1250.0_r8))*6.6e-5_r8/hcool
1554!
1555! Total cooling at the interface.
1556!
1557 qcool=lrad(i,j)+hsb+hlb-srad(i,j)*fc
1558 qbouy=tcff(i)*qcool+scff(i)*hlb*blk_cpw/hlv(i,j)
1559!
1560! Compute temperature and moisture change.
1561!
1562 IF ((qcool.gt.0.0_r8).and.(qbouy.gt.0.0_r8)) THEN
1563 lambd=6.0_r8/ &
1564 & (1.0_r8+(clam*qbouy/(wstar(i)+eps)**4)**0.75_r8)**r3
1565 hcool=lambd*blk_visw/(sqrt(rhoair(i)/rhosea(i))* &
1566 & wstar(i)+eps)
1567 deltc(i)=qcool*hcool/blk_tcw
1568 ELSE
1569 deltc(i)=0.0_r8
1570 END IF
1571 delqc(i)=cwet(i)*deltc(i)
1572# endif
1573 END DO
1574 END DO
1575!
1576
1577 DO i=istr-1,iendr
1578!
1579! Now complete the computation of the appropriate basic state quantities
1580! for this iteration.
1581!
1582# ifdef BBL_MODEL_NOT_YET
1583!
1584! Use wave info if we have it, two different options.
1585!
1586# ifdef WIND_WAVES
1587 zow(i)=(25._r8/pi)*lwave(i)*(wstar(i)/cwave(i))**4.5+ &
1588 & 0.11_r8*visair(i)/(wstar(i)+eps)
1589# else
1590 zow(i)=1200._r8*awave(i,j)*(awave(i,j)/lwave(i))**4.5+ &
1591 & 0.11_r8*visair(i)/(wstar(i)+eps)
1592# endif
1593# else
1594 zow(i)=charn(i)*wstar(i)*wstar(i)/g+ &
1595 & 0.11_r8*visair(i)/(wstar(i)+eps)
1596# endif
1597 rr(i)=zow(i)*wstar(i)/visair(i)
1598!
1599! Compute Monin-Obukhov stability parameter, Z/L.
1600!
1601 zoq(i)=min(1.15e-4_r8,5.5e-5_r8/rr(i)**0.6)
1602 zot(i)=zoq(i)
1603 zol(i)=vonkar*g*blk_zw(ng)* &
1604 & (tstar(i)*(1.0_r8+0.61_r8*q(i))+ &
1605 & 0.61_r8*tairk(i)*qstar(i))/ &
1606 & (tairk(i)*wstar(i)*wstar(i)* &
1607 & (1.0_r8+0.61_r8*q(i))+eps)
1608 l(i)=blk_zw(ng)/(zol(i)+eps)
1609!
1610! Evaluate stability functions at Z/L.
1611!
1612 wpsi(i)=bulk_psiu(zol(i),pi)
1613 tpsi(i)=bulk_psit(blk_zt(ng)/l(i),pi)
1614 qpsi(i)=bulk_psit(blk_zq(ng)/l(i),pi)
1615# ifdef COOL_SKIN
1616 delqc(i)=cwet(i)*deltc(i)
1617# endif
1618!
1619! Compute wind scaling parameters, Wstar.
1620!
1621 wstar1(i)=max(eps,delw(i)*vonkar/ &
1622 & (log(blk_zw(ng)/zow(i))-wpsi(i)))
1623 tstar1(i)=-(delt(i)-deltc(i))*vonkar/ &
1624 & (log(blk_zt(ng)/zot(i))-tpsi(i))
1625 qstar1(i)=-(delq(i)-delqc(i))*vonkar/ &
1626 & (log(blk_zq(ng)/zoq(i))-qpsi(i))
1627!
1628! Compute gustiness in wind speed.
1629!
1630 bf=-g/tairk(i)* &
1631 & wstar1(i)*(tstar1(i)+0.61_r8*tairk(i)*qstar1(i))
1632 IF (bf.gt.0.0_r8) THEN
1633 wgus(i)=blk_beta*(bf*blk_zabl)**r3
1634 ELSE
1635 wgus(i)=0.2_r8
1636 END IF
1637 delw1(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
1638# ifdef COOL_SKIN
1639!
1640!-----------------------------------------------------------------------
1641! Cool Skin correction.
1642!-----------------------------------------------------------------------
1643!
1644! Cool skin correction constants. Clam: part of Saunders constant
1645! lambda; Cwet: slope of saturation vapor.
1646!
1647 clam=16.0_r8*g*blk_cpw*(rhosea(i)*blk_visw)**3/ &
1648 & (blk_tcw*blk_tcw*rhoair(i)*rhoair(i))
1649!
1650! Set initial guesses for cool-skin layer thickness (Hcool).
1651!
1652 hcool=0.001_r8
1653!
1654! Backgound sensible and latent heat.
1655!
1656 hsb=-rhoair(i)*blk_cpa*wstar1(i)*tstar1(i)
1657 hlb=-rhoair(i)*hlv(i,j)*wstar1(i)*qstar1(i)
1658!
1659! Mean absoption in cool-skin layer.
1660!
1661 fc=0.065_r8+11.0_r8*hcool- &
1662 & (1.0_r8-exp(-hcool*1250.0_r8))*6.6e-5_r8/hcool
1663!
1664! Total cooling at the interface.
1665!
1666 qcool=lrad(i,j)+hsb+hlb-srad(i,j)*fc
1667 qbouy=tcff(i)*qcool+scff(i)*hlb*blk_cpw/hlv(i,j)
1668!
1669! Compute temperature and moisture change.
1670!
1671 IF ((qcool.gt.0.0_r8).and.(qbouy.gt.0.0_r8)) THEN
1672 lambd=6.0_r8/ &
1673 & (1.0_r8+(clam*qbouy/(wstar1(i)+eps)**4)**0.75_r8)**r3
1674 hcool=lambd*blk_visw/(sqrt(rhoair(i)/rhosea(i))* &
1675 & wstar1(i)+eps)
1676 deltc1(i)=qcool*hcool/blk_tcw
1677 ELSE
1678 deltc1(i)=0.0_r8
1679 END IF
1680 delqc(i)=cwet(i)*deltc1(i)
1681!
1682!-----------------------------------------------------------------------
1683! Adjoint Cool Skin correction.
1684!-----------------------------------------------------------------------
1685!
1686! Adjoint temperature and moisture change.
1687!
1688
1689
1690!^ tl_delQc(i)=tl_Cwet(i)*delTc1(i)+Cwet(i)*tl_delTc(i)
1691 ad_cwet(i)=ad_cwet(i)+deltc1(i)*ad_delqc(i)
1692 ad_deltc(i)=ad_deltc(i)+cwet(i)*ad_delqc(i)
1693 ad_delqc(i)=0.0_-r8
1694
1695 IF ((qcool.gt.0.0_r8).and.(qbouy.gt.0.0_r8)) THEN
1696 fac1=sqrt(rhoair(i)/rhosea(i))
1697 fac2=fac1*wstar1(i)+eps
1698 hcool=lambd*blk_visw/fac2
1699
1700!^ tl_delTc(i)=(tl_Qcool*Hcool+Qcool*tl_Hcool)/blk_tcw
1701 adfac=ad_deltc(i)/blk_tcw
1702 ad_qcool=ad_qcool+hcool*adfac
1703 ad_hcool=ad_hcool+qcool*adfac
1704 ad_deltc(i)=0.0_r8
1705
1706!^ tl_Hcool=tl_lambd*blk_visw/fac2-tl_fac2*Hcool/fac2
1707 adfac=ad_hcool/fac2
1708 ad_lambd=ad_lambd+blk_visw*adfac
1709 ad_fac2=-hcool*adfac
1710 ad_hcool=0.0_r8
1711
1712 IF (fac1.ne.0.0_r8) THEN
1713!^ tl_fac1=-0.5_r8*tl_rhoSea(i)*fac1/rhoSea(i)
1714 ad_rhosea(i)=ad_rhosea(i)-0.5*fac1*ad_fac1/rhosea(i)
1715 ad_fac1=0.0_r8
1716 ELSE
1717!^ tl_fac1=0.0_r8
1718 ad_fac1=0.0_r8
1719 END IF
1720
1721 fac1=(wstar1(i)+eps)**4
1722 fac2=clam*qbouy
1723 fac3=(fac2/fac1)**0.75_r8
1724 lambd=6.0_r8/(1.0_r8+fac3)**r3
1725
1726!^ tl_lambd=-r3*6.0_r8*tl_fac3/(1.0_r8+fac3)**(r3+1.0_r8)
1727 ad_fac3=-r3*6.0_r8*ad_lambd/(1.0_r8+fac3)**(r3+1.0_r8)
1728 ad_lambd=0.0_r8
1729
1730!^ tl_fac3=0.75_r8*(fac2/fac1)**(-0.25_r8)* &
1731!^ & (tl_fac2/fac1-tl_fac1*fac2/(fac1*fac1))
1732 adfac=0.75_r8*(fac2/fac1)**(-0.25_r8)*ad_fac3
1733 ad_fac2=adfac/fac1
1734 ad_fac1=-fac2*adfac/(fac1*fac1)
1735 ad_fac3=0.0_r8
1736
1737!^ tl_fac2=tl_Clam*Qbouy+Clam*tl_Qbouy
1738
1739 ad_clam=ad_clam+qbouy*ad_fac2
1740 ad_qbouy=ad_qbouy+clam*ad_fac2
1741 ad_fac2=0.0_r8
1742
1743!^ tl_fac1=4.0_r8*tl_Wstar(i)*(Wstar1(i)+eps)**3
1744
1745 ad_wstar(i)=ad_wstar(i)+4.0_r8*ad_fac1*(wstar1(i)+eps)**3
1746 ad_fac1=0.0_r8
1747
1748 ELSE
1749!^ tl_delTc(i)=0.0_r8
1750 ad_deltc(i)=0.0_r8
1751 END IF
1752!
1753! Adjoint of Total cooling at the interface.
1754!
1755 clam=16.0_r8*g*blk_cpw*(rhosea(i)*blk_visw)**3/ &
1756 & (blk_tcw*blk_tcw*rhoair(i)*rhoair(i))
1757!
1758! Set initial guesses for cool-skin layer thickness (Hcool).
1759!
1760 hcool=0.001_r8
1761!
1762! Backgound sensible and latent heat.
1763!
1764 hsb=-rhoair(i)*blk_cpa*wstar1(i)*tstar1(i)
1765 hlb=-rhoair(i)*hlv(i,j)*wstar1(i)*qstar1(i)
1766!
1767! Mean absoption in cool-skin layer.
1768!
1769 fc=0.065_r8+11.0_r8*hcool- &
1770 & (1.0_r8-exp(-hcool*1250.0_r8))*6.6e-5_r8/hcool
1771!
1772! Total cooling at the interface.
1773!
1774 qcool=lrad(i,j)+hsb+hlb-srad(i,j)*fc
1775 qbouy=tcff(i)*qcool+scff(i)*hlb*blk_cpw/hlv(i,j)
1776
1777!^ tl_Qbouy=tl_Tcff(i)*Qcool+Tcff(i)*tl_Qcool+ &
1778!^ & (tl_Scff(i)*Hlb*blk_Cpw+ &
1779!^ & Scff(i)*tl_Hlb*blk_Cpw)/Hlv(i,j)- &
1780!^ & tl_Hlv(i,j)*Scff(i)*Hlb*blk_Cpw/(Hlv(i,j)*Hlv(i,j))
1781
1782 adfac=ad_qbouy*blk_cpw/hlv(i,j)
1783 ad_tcff(i)=ad_tcff(i)+qcool*ad_qbouy
1784 ad_qcool=ad_qcool+tcff(i)*ad_qbouy
1785 ad_scff(i)=ad_scff(i)+hlb*adfac
1786 ad_hlb=ad_hlb+scff(i)*adfac
1787 ad_hlv(i,j)=ad_hlv(i,j)- &
1788 & scff(i)*hlb*adfac/hlv(i,j)
1789 ad_qbouy=0.0_r8
1790
1791
1792!^ tl_Qcool=tl_LRad(i,j)+tl_Hsb+tl_Hlb
1793 ad_lrad(i,j)=ad_lrad(i,j)+ad_qcool
1794 ad_hsb=ad_hsb+ad_qcool
1795 ad_hlb=ad_hlb+ad_qcool
1796 ad_qcool=0.0_r8
1797!
1798! Adjoint Backgound sensible and latent heat.
1799!
1800!^ tl_Hlb=-rhoAir(i)*(tl_Hlv(i,j)*Wstar1(i)*Qstar1(i)+ &
1801!^ & Hlv(i,j)*(tl_Wstar(i)*Qstar1(i)+ &
1802!^ & Wstar1(i)*tl_Qstar(i))
1803
1804 adfac=-rhoair(i)*ad_hlb
1805 adfac1=hlv(i,j)*adfac
1806 ad_hlv(i,j)=ad_hlv(i,j)+wstar1(i)*qstar1(i)*adfac
1807 ad_wstar(i)=ad_wstar(i)+qstar1(i)*adfac1
1808 ad_qstar(i)=ad_qstar(i)+wstar1(i)*adfac1
1809 ad_hlb=0.0_r8
1810
1811!^ tl_Hsb=-rhoAir(i)*blk_Cpa*(tl_Wstar(i)*Tstar1(i)+ &
1812!^ & Wstar1(i)*tl_Tstar(i))
1813 adfac=-rhoair(i)*blk_cpa*ad_hsb
1814 ad_wstar(i)=ad_wstar(i)+tstar1(i)*adfac
1815 ad_tstar(i)=ad_tstar(i)+wstar1(i)*adfac
1816 ad_hsb=0.0_r8
1817
1818!
1819! Adjoint Cool skin correction constants. Clam: part of Saunders constant
1820! lambda; Cwet: slope of saturation vapor.
1821!
1822!^ tl_Clam=48.0_r8*g*blk_Cpw*tl_rhoSea(i)*blk_visw* &
1823!^ & (rhoSea(i)*blk_visw)**2/ &
1824!^ & (blk_tcw*blk_tcw*rhoAir(i)*rhoAir(i))
1825
1826 ad_rhosea(i)=ad_rhosea(i)+ &
1827 & 48.0_r8*g*blk_cpw*ad_clam*blk_visw* &
1828 & (rhosea(i)*blk_visw)**2/ &
1829 & (blk_tcw*blk_tcw*rhoair(i)*rhoair(i))
1830 ad_clam=0.0_r8
1831# endif
1832
1833!
1834! Adjoint Compute gustiness in wind speed.
1835!
1836 IF (delw1(i).ne.0.0_r8)THEN
1837!^ tl_delW(i)=tl_Wgus(i)*Wgus(i)/delW1(i)
1838 ad_wgus(i)=ad_wgus(i)+wgus(i)*ad_delw(i)/delw1(i)
1839 ad_delw(i)=0.0_r8
1840 ELSE
1841!^ tl_delW(i)=0.0_r8
1842 ad_delw(i)=0.0_r8
1843 END IF
1844
1845 IF (bf.gt.0.0_r8) THEN
1846!^ tl_Wgus(i)=r3*blk_beta*tl_Bf*blk_Zabl* &
1847!^ & (Bf*blk_Zabl)**(r3-1.0_r8)
1848 ad_bf=ad_bf+r3*blk_beta*ad_wgus(i)*blk_zabl* &
1849 & (bf*blk_zabl)**(r3-1.0_r8)
1850 ad_wgus(i)=0.0_r8
1851 ELSE
1852!^ tl_Wgus(i)=0.0_r8
1853 ad_wgus(i)=0.0_r8
1854 END IF
1855
1856!^ tl_Bf=-g/TairK(i)*( &
1857!^ & tl_Wstar(i)*(Tstar1(i)+0.61_r8*TairK(i)*Qstar1(i))+ &
1858!^ & Wstar1(i)*(tl_Tstar(i)+0.61_r8*TairK(i)*tl_Qstar(i)) )
1859 adfac=-g/tairk(i)*ad_bf
1860 adfac1=wstar1(i)*adfac
1861 ad_wstar(i)=ad_wstar(i)+ &
1862 & (tstar1(i)+0.61_r8*tairk(i)*qstar1(i))*adfac
1863 ad_tstar(i)=ad_tstar(i)+adfac1
1864 ad_qstar(i)=ad_qstar(i)+0.61_r8*tairk(i)*adfac1
1865 ad_bf=0.0_r8
1866
1867!
1868! Adjoint of wind scaling parameters, Wstar.
1869!
1870
1871 fac1=blk_zq(ng)/zoq(i)
1872 fac2=log(fac1)
1873
1874!^ tl_Qstar(i)=-(tl_delQ(i)-tl_delQc(i))*vonKar/(fac2-Qpsi(i))-&
1875!^ & (tl_fac2-tl_Qpsi(i))*Qstar1(i)/(fac2-Qpsi(i))
1876 adfac=-ad_qstar(i)*vonkar/(fac2-qpsi(i))
1877 adfac1=-ad_qstar(i)*qstar1(i)/(fac2-qpsi(i))
1878 ad_delq(i)=ad_delq(i)+adfac
1879 ad_delqc(i)=ad_delqc(i)-adfac
1880 ad_fac2=adfac1
1881 ad_qpsi(i)=ad_qpsi(i)-adfac1
1882 ad_qstar(i)=0.0_r8
1883
1884!^ tl_fac2=tl_fac1/fac1
1885 ad_fac1=ad_fac2/fac1
1886 ad_fac2=0.0_r8
1887
1888!^ tl_fac1=-tl_ZoQ(i)*fac1/ZoQ(i)
1889 ad_zoq(i)=ad_zoq(i)-fac1*ad_fac1/zoq(i)
1890 ad_fac1=0.0_r8
1891
1892 fac1=blk_zt(ng)/zot(i)
1893 fac2=log(fac1)
1894
1895!^ tl_Tstar(i)=-(tl_delT(i)-tl_delTc(i))*vonKar/(fac2-Tpsi(i))-&
1896!^ & (tl_fac2-tl_Tpsi(i))*Tstar1(i)/(fac2-Tpsi(i))
1897 adfac=-ad_tstar(i)*vonkar/(fac2-tpsi(i))
1898 adfac1=-ad_tstar(i)*tstar1(i)/(fac2-tpsi(i))
1899 ad_delt(i)=ad_delt(i)+adfac
1900 ad_deltc(i)=ad_deltc(i)-adfac
1901 ad_fac2=adfac1
1902 ad_tpsi(i)=ad_tpsi(i)-adfac1
1903 ad_tstar(i)=0.0_r8
1904
1905!^ tl_fac2=tl_fac1/fac1
1906 ad_fac1=ad_fac1+ad_fac2/fac1
1907 ad_fac2=0.0_r8
1908
1909!^ tl_fac1=-tl_ZoT(i)*fac1/ZoT(i)
1910 ad_zot(i)=ad_zot(i)-fac1*ad_fac1/zot(i)
1911 ad_fac1=0.0_r8
1912
1913 fac1=blk_zw(ng)/zow(i)
1914 fac2=log(fac1)
1915 fac3=delw(i)*vonkar/(fac2-wpsi(i))
1916
1917!^ tl_Wstar(i)=(0.5_r8-SIGN(0.5_r8,eps-fac3))*tl_fac3
1918 ad_fac3=ad_wstar(i)*(0.5_r8-sign(0.5_r8,eps-fac3))
1919 ad_wstar(i)=0.0_r8
1920
1921!^ tl_fac3=tl_delW(i)*vonKar/(fac2-Wpsi(i))- &
1922!^ & (tl_fac2-tl_Wpsi(i))*fac3/(fac2-Wpsi(i))
1923 adfac=ad_fac3/(fac2-wpsi(i))
1924 adfac1=fac3*adfac
1925 ad_delw(i)=ad_delw(i)+vonkar*adfac
1926 ad_fac2=-adfac1
1927 ad_wpsi(i)=ad_wpsi(i)+adfac1
1928 ad_fac3=0.0_r8
1929
1930!^ tl_fac2=tl_fac1/fac1
1931 ad_fac1=ad_fac1+ad_fac2/fac1
1932 ad_fac2=0.0_r8
1933
1934!^ tl_fac1=-tl_ZoW(i)*fac1/ZoW(i)
1935 ad_zow(i)=ad_zow(i)-fac1*ad_fac1/zow(i)
1936 ad_fac1=0.0_r8
1937
1938# ifdef COOL_SKIN
1939!^ tl_delQc(i)=tl_Cwet(i)*delTc(i)+Cwet(i)*tl_delTc(i)
1940 ad_cwet(i)=ad_cwet(i)+deltc(i)*ad_delqc(i)
1941 ad_deltc(i)=ad_deltc(i)+cwet(i)*ad_delqc(i)
1942 ad_delqc(i)=0.0_r8
1943
1944!^ tl_Cwet(i)=0.622_r8*(tl_Hlv(i,j)*Qsea(i)+ &
1945!^ & Hlv(i,j)*tl_Qsea(i))/ &
1946!^ & (blk_Rgas*TseaK(i)*TseaK(i))- &
1947!^ & 2.0_r8*blk_Rgas*tl_TseaK(i)*TseaK(i)*Cwet(i)/ &
1948!^ & (blk_Rgas*TseaK(i)*TseaK(i))
1949
1950 adfac=ad_cwet(i)/(blk_rgas*tseak(i)*tseak(i))
1951 adfac1=2.0_r8*blk_rgas*tseak(i)*cwet(i)*adfac
1952 adfac2=0.622_r8*adfac
1953 ad_hlv(i,j)=ad_hlv(i,j)+qsea(i)*adfac2
1954 ad_qsea(i)=ad_qsea(i)+hlv(i,j)*adfac2
1955 ad_tseak(i)=ad_tseak(i)-adfac1
1956 ad_cwet(i)=0.0_r8
1957
1958# endif
1959
1960!
1961! Adjoint stability functions at Z/L.
1962!
1963 fac=blk_zq(ng)/l(i)
1964!^ tl_Qpsi(i)=tl_bulk_psit(tl_fac,fac,pi)
1965 ad_fac=ad_bulk_psit(ad_qpsi(i),fac,pi)
1966 ad_qpsi(i)=0.0_r8
1967
1968!^ tl_fac=-tl_L(i)*fac/L(i)
1969 ad_l(i)=ad_l(i)-ad_fac*fac/l(i)
1970 ad_fac=0.0_r8
1971
1972 fac=blk_zt(ng)/l(i)
1973!^ tl_Tpsi(i)=tl_bulk_psit(tl_fac,fac,pi)
1974 ad_fac=ad_bulk_psit(ad_tpsi(i),fac,pi)
1975 ad_tpsi(i)=0.0_r8
1976
1977!^ tl_fac=-tl_L(i)*fac/L(i)
1978 ad_l(i)=ad_l(i)-ad_fac*fac/l(i)
1979 ad_fac=0.0_r8
1980
1981!^ tl_Wpsi(i)=tl_bulk_psiu(tl_ZoL(i),ZoL(i),pi)
1982 ad_zol(i)=ad_zol(i)+ad_bulk_psiu(ad_wpsi(i),zol(i),pi)
1983 ad_wpsi(i)=0.0_r8
1984
1985!
1986! Adjoint Monin-Obukhov stability parameter, Z/L.
1987!
1988!^ tl_L(i)=-tl_ZoL(i)*L(i)/(ZoL(i)+eps)
1989 ad_zol(i)=ad_zol(i)-l(i)*ad_l(i)/(zol(i)+eps)
1990 ad_l(i)=0.0_r8
1991
1992!^ tl_ZoL(i)=vonKar*g*blk_ZW(ng)* &
1993!^ & (tl_Tstar(i)*(1.0_r8+0.61_r8*Q(i))+ &
1994!^ & 0.61_r8*TairK(i)*tl_Qstar(i))/ &
1995!^ & (TairK(i)*Wstar(i)*Wstar(i)* &
1996!^ & (1.0_r8+0.61_r8*Q(i))+eps)- &
1997!^ & 2.0_r8*TairK(i)*tl_Wstar(i)*Wstar(i)* &
1998!^ & (1.0_r8+0.61_r8*Q(i))*ZoL(i)/ &
1999!^ & (TairK(i)*Wstar(i)*Wstar(i)* &
2000!^ & (1.0_r8+0.61_r8*Q(i))+eps) )
2001 adfac=vonkar*g*blk_zw(ng)*ad_zol(i)/ &
2002 & (tairk(i)*wstar(i)*wstar(i)* &
2003 & (1.0_r8+0.61_r8*q(i))+eps)
2004 ad_tstar(i)=ad_tstar(i)+(1.0_r8+0.61_r8*q(i))*adfac
2005 ad_qstar(i)=ad_qstar(i)+0.61_r8*tairk(i)*adfac
2006 ad_wstar(i)=ad_wstar(i)-2.0_r8*tairk(i)*wstar(i)* &
2007 & (1.0_r8+0.61_r8*q(i))*zol(i)*ad_zol(i)/ &
2008 & (tairk(i)*wstar(i)*wstar(i)* &
2009 & (1.0_r8+0.61_r8*q(i))+eps)
2010 ad_zol(i)=0.0_r8
2011
2012!^ tl_ZoT(i)=tl_ZoQ(i)
2013 ad_zoq(i)=ad_zoq(i)+ad_zot(i)
2014 ad_zot(i)=0.0_r8
2015
2016!^ tl_ZoQ(i)= &
2017!^ & -(0.5_r8-SIGN(0.5_r8,5.5e-5_r8/Rr(i)**0.6-1.15e-4_r8)) &
2018!^ & *0.6_r8*5.5e-5_r8*tl_Rr(i)/Rr(i)**1.6
2019 ad_rr(i)=ad_rr(i)- &
2020 & (0.5_r8-sign(0.5_r8,5.5e-5_r8/rr(i)**0.6-1.15e-4_r8)) &
2021 & *0.6_r8*5.5e-5_r8*ad_zoq(i)/rr(i)**1.6
2022 ad_zoq(i)=0.0_r8
2023
2024!^ tl_Rr(i)=(tl_ZoW(i)*Wstar(i)+ZoW(i)*tl_Wstar(i))/VisAir(i)
2025 adfac=ad_rr(i)/visair(i)
2026 ad_zow(i)=ad_zow(i)+wstar(i)*adfac
2027 ad_wstar(i)=ad_wstar(i)+zow(i)*adfac
2028 ad_rr(i)=0.0_r8
2029
2030# ifdef BBL_MODEL_NOT_YET
2031# ifdef WIND_WAVES
2032!^ tl_ZoW(i)=(25._r8/pi)*Lwave(i)*4.5_r8*tl_Wstar(i)* &
2033!^ & (Wstar(i)/Cwave(i))**3.5/Cwave(i)- &
2034!^ & tl_Wstar(i)*0.11_r8*VisAir(i)/ &
2035!^ & ((Wstar(i)+eps)*(Wstar(i)+eps))
2036 ad_wstar(i)=ad_wstar(i)+(25._r8/pi)*lwave(i)*4.5_r8* &
2037 & (wstar(i)/cwave(i))**3.5*ad_zow(i)/cwave(i)- &
2038 & 0.11_r8*visair(i)*ad_zow(i)/ &
2039 & ((wstar(i)+eps)*(wstar(i)+eps))
2040 ad_zow(i)0.0_r8
2041# else
2042!^ tl_ZoW(i)=-tl_Wstar(i)*0.11_r8*VisAir(i)/ &
2043!^ & ((Wstar(i)+eps)*(Wstar(i)+eps))
2044 ad_wstar(i)=ad_wstar(i)-0.11_r8*visair(i)*ad_zow(i)/ &
2045 & ((wstar(i)+eps)*(wstar(i)+eps))
2046 ad_zow(i)0.0_r8
2047# endif
2048# else
2049!^ tl_ZoW(i)=2.0_r8*charn(i)*tl_Wstar(i)*Wstar(i)/g- &
2050!^ & tl_Wstar(i)*0.11_r8*VisAir(i)/ &
2051!^ & ((Wstar(i)+eps)*(Wstar(i)+eps))
2052 adfac=ad_zow(i)/g
2053 ad_wstar(i)=ad_wstar(i)+2.0_r8*charn(i)*wstar(i)*adfac- &
2054 & 0.11_r8*visair(i)*ad_zow(i)/ &
2055 & ((wstar(i)+eps)*(wstar(i)+eps))
2056 ad_zow(i)=0.0_r8
2057
2058# endif
2059 END DO
2060 END DO
2061!
2062! End of Adjoint Iteration Loop
2063!
2064
2065 DO i=istr-1,iendr
2066!
2067! Input bulk parameterization fields.
2068!
2069 wmag(i)=sqrt(uwind(i,j)*uwind(i,j)+vwind(i,j)*vwind(i,j))
2070 pairm=pair(i,j)
2071 tairc(i)=tair(i,j)
2072 tairk(i)=tairc(i)+273.16_r8
2073 tseac(i)=t(i,j,n(ng),nrhs,itemp)
2074 tseak(i)=tseac(i)+273.16_r8
2075 rhosea(i)=rho(i,j,n(ng))+1000.0_r8
2076 rh=hair(i,j)
2077 srad(i,j)=srflx(i,j)*hscale
2078 tcff(i)=alpha(i,j)
2079 scff(i)=beta(i,j)
2080!
2081! Initialize.
2082!
2083 deltc(i)=0.0_r8
2084 delqc(i)=0.0_r8
2085 lheat(i,j)=lhflx(i,j)*hscale
2086 sheat(i,j)=shflx(i,j)*hscale
2087 taur=0.0_r8
2088 taux(i,j)=0.0_r8
2089 tauy(i,j)=0.0_r8
2090!
2091!-----------------------------------------------------------------------
2092! Compute net longwave radiation (W/m2), LRad.
2093!-----------------------------------------------------------------------
2094
2095# if defined LONGWAVE
2096!
2097! Use Berliand (1952) formula to calculate net longwave radiation.
2098! The equation for saturation vapor pressure is from Gill (Atmosphere-
2099! Ocean Dynamics, pp 606). Here the coefficient in the cloud term
2100! is assumed constant, but it is a function of latitude varying from
2101! 1.0 at poles to 0.5 at the equator).
2102!
2103 cff=(0.7859_r8+0.03477_r8*tairc(i))/ &
2104 & (1.0_r8+0.00412_r8*tairc(i))
2105 e_sat=10.0_r8**cff ! saturation vapor pressure (hPa or mbar)
2106 vap_p=e_sat*rh ! water vapor pressure (hPa or mbar)
2107 cff2=tairk(i)*tairk(i)*tairk(i)
2108 cff1=cff2*tairk(i)
2109 lrad(i,j)=-emmiss*stefbo* &
2110 & (cff1*(0.39_r8-0.05_r8*sqrt(vap_p))* &
2111 & (1.0_r8-0.6823_r8*cloud(i,j)*cloud(i,j))+ &
2112 & cff2*4.0_r8*(tseak(i)-tairk(i)))
2113
2114# elif defined LONGWAVE_OUT
2115!
2116! Treat input longwave data as downwelling radiation only and add
2117! outgoing IR from model sea surface temperature.
2118!
2119 lrad(i,j)=lrflx(i,j)*hscale- &
2120 & emmiss*stefbo*tseak(i)*tseak(i)*tseak(i)*tseak(i)
2121
2122# else
2123 lrad(i,j)=lrflx(i,j)*hscale
2124# endif
2125# ifdef MASKING
2126 lrad(i,j)=lrad(i,j)*rmask(i,j)
2127# endif
2128!
2129!-----------------------------------------------------------------------
2130! Compute specific humidities (kg/kg).
2131!
2132! note that Qair(i) is the saturation specific humidity at Tair
2133! Q(i) is the actual specific humidity
2134! Qsea(i) is the saturation specific humidity at Tsea
2135!
2136! Saturation vapor pressure in mb is first computed and then
2137! converted to specific humidity in kg/kg
2138!
2139! The saturation vapor pressure is computed from Teten formula
2140! using the approach of Buck (1981):
2141!
2142! Esat(mb) = (1.0007_r8+3.46E-6_r8*PairM(mb))*6.1121_r8*
2143! EXP(17.502_r8*TairC(C)/(240.97_r8+TairC(C)))
2144!
2145! The ambient vapor is found from the definition of the
2146! Relative humidity:
2147!
2148! RH = W/Ws*100 ~ E/Esat*100 E = RH/100*Esat if RH is in %
2149! E = RH*Esat if RH fractional
2150!
2151! The specific humidity is then found using the relationship:
2152!
2153! Q = 0.622 E/(P + (0.622-1)e)
2154!
2155! Q(kg/kg) = 0.62197_r8*(E(mb)/(PairM(mb)-0.378_r8*E(mb)))
2156!
2157!-----------------------------------------------------------------------
2158!
2159! Compute air saturation vapor pressure (mb), using Teten formula.
2160!
2161 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
2162 & exp(17.502_r8*tairc(i)/(240.97_r8+tairc(i)))
2163!
2164! Compute specific humidity at Saturation, Qair (kg/kg).
2165!
2166 qair(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff))
2167!
2168! Compute specific humidity, Q (kg/kg).
2169!
2170 IF (rh.lt.2.0_r8) THEN !RH fraction
2171 cff=cff*rh !Vapor pres (mb)
2172 q(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff)) !Spec hum (kg/kg)
2173 ELSE !RH input was actually specific humidity in g/kg
2174 q(i)=rh/1000.0_r8 !Spec Hum (kg/kg)
2175 END IF
2176!
2177! Compute water saturation vapor pressure (mb), using Teten formula.
2178!
2179 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
2180 & exp(17.502_r8*tseac(i)/(240.97_r8+tseac(i)))
2181!
2182! Vapor Pressure reduced for salinity (Kraus and Businger, 1994, pp42).
2183!
2184 cff=cff*0.98_r8
2185!
2186! Compute Qsea (kg/kg) from vapor pressure.
2187!
2188 qsea(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff))
2189!
2190!-----------------------------------------------------------------------
2191! Compute Monin-Obukhov similarity parameters for wind (Wstar),
2192! heat (Tstar), and moisture (Qstar), Liu et al. (1979).
2193!-----------------------------------------------------------------------
2194!
2195! Moist air density (kg/m3).
2196!
2197 rhoair(i)=pairm*100.0_r8/(blk_rgas*tairk(i)* &
2198 & (1.0_r8+0.61_r8*q(i)))
2199!
2200! Kinematic viscosity of dry air (m2/s), Andreas (1989).
2201!
2202 visair(i)=1.326e-5_r8* &
2203 & (1.0_r8+tairc(i)*(6.542e-3_r8+tairc(i)* &
2204 & (8.301e-6_r8-4.84e-9_r8*tairc(i))))
2205!
2206! Compute latent heat of vaporization (J/kg) at sea surface, Hlv.
2207!
2208 hlv(i,j)=(2.501_r8-0.00237_r8*tseac(i))*1.0e+6_r8
2209!
2210! Assume that wind is measured relative to sea surface and include
2211! gustiness.
2212!
2213 wgus(i)=0.5_r8
2214 delw(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
2215 delq(i)=qsea(i)-q(i)
2216 delt(i)=tseac(i)-tairc(i)
2217!
2218! Neutral coefficients.
2219!
2220 zow(i)=0.0001_r8
2221 u10(i)=delw(i)*log(10.0_r8/zow(i))/log(blk_zw(ng)/zow(i))
2222 wstar(i)=0.035_r8*u10(i)
2223 zo10(i)=0.011_r8*wstar(i)*wstar(i)/g+ &
2224 & 0.11_r8*visair(i)/wstar(i)
2225 cd10(i)=(vonkar/log(10.0_r8/zo10(i)))**2
2226 ch10(i)=0.00115_r8
2227 ct10(i)=ch10(i)/sqrt(cd10(i))
2228 zot10(i)=10.0_r8/exp(vonkar/ct10(i))
2229 cd=(vonkar/log(blk_zw(ng)/zo10(i)))**2
2230!
2231! Compute Richardson number.
2232!
2233 ct(i)=vonkar/log(blk_zt(ng)/zot10(i)) ! T transfer coefficient
2234 cc(i)=vonkar*ct(i)/cd
2235 deltc(i)=0.0_r8
2236# ifdef COOL_SKIN
2237 deltc(i)=blk_dter
2238# endif
2239 ribcu(i)=-blk_zw(ng)/(blk_zabl*0.004_r8*blk_beta**3)
2240 ri(i)=-g*blk_zw(ng)*((delt(i)-deltc(i))+ &
2241 & 0.61_r8*tairk(i)*delq(i))/ &
2242 & (tairk(i)*delw(i)*delw(i))
2243 IF (ri(i).lt.0.0_r8) THEN
2244 zetu(i)=cc(i)*ri(i)/(1.0_r8+ri(i)/ribcu(i)) ! Unstable
2245 ELSE
2246 zetu(i)=cc(i)*ri(i)/(1.0_r8+3.0_r8*ri(i)/cc(i)) ! Stable
2247 END IF
2248 l10(i)=blk_zw(ng)/zetu(i)
2249!
2250! First guesses for Monon-Obukhov similarity scales.
2251!
2252 wstar(i)=delw(i)*vonkar/(log(blk_zw(ng)/zo10(i))- &
2253 & bulk_psiu(blk_zw(ng)/l10(i),pi))
2254 tstar(i)=-(delt(i)-deltc(i))*vonkar/ &
2255 & (log(blk_zt(ng)/zot10(i))- &
2256 & bulk_psit(blk_zt(ng)/l10(i),pi))
2257 qstar(i)=-(delq(i)-delqc(i))*vonkar/ &
2258 & (log(blk_zq(ng)/zot10(i))- &
2259 & bulk_psit(blk_zq(ng)/l10(i),pi))
2260!
2261! Modify Charnock for high wind speeds. The 0.125 factor below is for
2262! 1.0/(18.0-10.0).
2263!
2264 IF (delw(i).gt.18.0_r8) THEN
2265 charn(i)=0.018_r8
2266 ELSE IF ((10.0_r8.lt.delw(i)).and.(delw(i).le.18.0_r8)) THEN
2267 charn(i)=0.011_r8+ &
2268 & 0.125_r8*(0.018_r8-0.011_r8)*(delw(i)-10._r8)
2269 ELSE
2270 charn(i)=0.011_r8
2271 END IF
2272
2273
2274!
2275! Adjoint of Charnock for high wind speeds. The 0.125 factor below is for
2276! 1.0/(18.0-10.0).
2277!
2278
2279 IF (delw(i).gt.18.0_r8) THEN
2280!^ tl_charn(i)=0.0_r8
2281 ad_charn(i)=0.0_r8
2282 ELSE IF ((10.0_r8.lt.delw(i)).and.(delw(i).le.18.0_r8)) THEN
2283!^ tl_charn(i)=0.0_r8
2284 ad_charn(i)=0.0_r8
2285 ELSE
2286!^ tl_charn(i)=0.0_r8
2287 ad_charn(i)=0.0_r8
2288 END IF
2289
2290!
2291! Adjoint First guesses for Monon-Obukhov similarity scales.
2292!
2293
2294 fac1=blk_zq(ng)/l10(i)
2295 fac2=log(blk_zq(ng)/zot10(i))
2296
2297!^ tl_Qstar(i)=-(tl_delQ(i)-tl_delQc(i))*vonKar/ &
2298!^ & (fac2-bulk_psit(fac1,pi))- &
2299!^ & (tl_fac2-tl_bulk_psit(tl_fac1,fac1,pi))*Qstar(i)/ &
2300!^ & (fac2-bulk_psit(fac1,pi))
2301
2302 adfac=ad_qstar(i)*vonkar/ &
2303 & (fac2-bulk_psit(fac1,pi))
2304 cff=qstar(i)/(fac2-bulk_psit(fac1,pi))
2305 ad_delq(i)=ad_delq(i)-adfac
2306 ad_delqc(i)=ad_delqc(i)+adfac
2307 ad_fac2=-cff*ad_qstar(i)
2308!^ ad_fac1=ad_bulk_psit(ad_Qstar(i),fac1,pi)*cff
2309 ad_fac1=ad_bulk_psit(cff*ad_qstar(i),fac1,pi)
2310 ad_qstar(i)=0.0_r8
2311
2312!^ tl_fac2=-tl_ZoT10(i)/ZoT10(i)
2313 ad_zot10(i)=ad_zot10(i)-ad_fac2/zot10(i)
2314 ad_fac2=0.0_r8
2315
2316!^ tl_fac1=-tl_L10(i)*fac1/L10(i)
2317 ad_l10(i)=ad_l10(i)-ad_fac1*fac1/l10(i)
2318 ad_fac1=0.0_r8
2319
2320 fac1=blk_zt(ng)/l10(i)
2321 fac2=log(blk_zt(ng)/zot10(i))
2322
2323!^ tl_Tstar(i)=-(tl_delT(i)-tl_delTc(i))*vonKar/ &
2324!^ & (fac2-bulk_psit(fac1,pi))- &
2325!^ & (tl_fac2-tl_bulk_psit(tl_fac1,fac1,pi))*Tstar(i)/ &
2326!^ & (fac2-bulk_psit(fac1,pi))
2327 adfac=ad_tstar(i)*vonkar/ &
2328 & (fac2-bulk_psit(fac1,pi))
2329 cff=tstar(i)/(fac2-bulk_psit(fac1,pi))
2330 ad_delt(i)=ad_delt(i)-adfac
2331 ad_deltc(i)=ad_deltc(i)+adfac
2332 ad_fac2=-cff*ad_tstar(i)
2333!^ ad_fac1=ad_bulk_psit(ad_Tstar(i),fac1,pi)*cff
2334 ad_fac1=ad_bulk_psit(cff*ad_tstar(i),fac1,pi)
2335 ad_tstar(i)=0.0_r8
2336
2337!^ tl_fac2=-tl_ZoT10(i)/ZoT10(i)
2338 ad_zot10(i)=ad_zot10(i)-ad_fac2/zot10(i)
2339 ad_fac2=0.0_r8
2340
2341!^ tl_fac1=-tl_L10(i)*fac1/L10(i)
2342 ad_l10(i)=ad_l10(i)-fac1*ad_fac1/l10(i)
2343 ad_fac1=0.0_r8
2344
2345 fac1=blk_zw(ng)/l10(i)
2346 fac2=log(blk_zw(ng)/zo10(i))
2347
2348!^ tl_Wstar(i)=-(tl_fac2-tl_bulk_psiu(tl_fac1,fac1,pi))*Wstar(i) &
2349!^ & /(fac2-bulk_psiu(fac1,pi))
2350 cff=wstar(i)/(fac2-bulk_psiu(fac1,pi))
2351 ad_fac2=-cff*ad_wstar(i)
2352!^ ad_fac1=cff*ad_bulk_psiu(ad_Wstar(i),fac1,pi)
2353 ad_fac1=ad_bulk_psiu(cff*ad_wstar(i),fac1,pi)
2354 ad_wstar(i)=0.0_r8
2355
2356!^ tl_fac2=-tl_Zo10(i)/Zo10(i)
2357 ad_zo10(i)=ad_zo10(i)-ad_fac2/zo10(i)
2358 ad_fac2=0.0_r8
2359
2360!^ tl_fac1=-tl_L10(i)*fac1/L10(i)
2361 ad_l10(i)=ad_l10(i)-ad_fac1*fac1/l10(i)
2362 ad_fac1=0.0_r8
2363
2364!
2365! Adjoint Richardson number.
2366!
2367
2368!^ tl_L10(i)=-L10(i)*L10(i)*tl_Zetu(i)/blk_ZW(ng)
2369 ad_zetu(i)=ad_zetu(i)-l10(i)*l10(i)*ad_l10(i)/blk_zw(ng)
2370 ad_l10(i)=0.0_r8
2371
2372 IF (ri(i).lt.0.0_r8) THEN
2373!^ tl_Zetu(i)=(tl_CC(i)*Ri(i)+CC(i)*tl_Ri(i))/ &
2374!^ & (1.0_r8+Ri(i)/Ribcu(i))- &
2375!^ & (tl_Ri(i)/Ribcu(i))*Zetu(i)/(1.0_r8+Ri(i)/Ribcu(i))
2376 adfac=ad_zetu(i)/(1.0_r8+ri(i)/ribcu(i))
2377 ad_cc(i)=ad_cc(i)+ri(i)*adfac
2378 ad_ri(i)=ad_ri(i)+cc(i)*adfac-zetu(i)*adfac/ribcu(i)
2379 ad_zetu(i)=0.0_r8
2380 ELSE
2381 fac=3.0_r8*ri(i)/cc(i)
2382!^ tl_Zetu(i)=(tl_CC(i)*Ri(i)+CC(i)*tl_Ri(i))/(1.0_r8+fac) &
2383!^ & -tl_fac*Zetu(i)/(1.0_r8+fac)
2384 adfac=ad_zetu(i)/(1.0_r8+fac)
2385 ad_cc(i)=ad_cc(i)+ri(i)*adfac
2386 ad_ri(i)=ad_ri(i)+cc(i)*adfac
2387 ad_fac=-zetu(i)*adfac
2388 ad_zetu(i)=0.0_r8
2389
2390!^ tl_fac=3.0_r8*tl_Ri(i)/CC(i)-tl_CC(i)*fac/CC(i)
2391 ad_ri(i)=ad_ri(i)+3.0_r8*ad_fac/cc(i)
2392 ad_cc(i)=ad_cc(i)-ad_fac*fac/cc(i)
2393 ad_fac=0.0_r8
2394 END IF
2395
2396 fac=1.0/(tairk(i)*delw(i)*delw(i))
2397
2398!^ tl_Ri(i)=-fac*g*blk_ZW(ng)*((tl_delT(i)-tl_delTc(i))+ &
2399!^ & 0.61_r8*TairK(i)*tl_delQ(i))
2400!^
2401 adfac=-fac*g*blk_zw(ng)*ad_ri(i)
2402 ad_delt(i)=ad_delt(i)+adfac
2403 ad_deltc(i)=ad_deltc(i)-adfac
2404 ad_delq(i)=ad_delq(i)+0.61_r8*tairk(i)*adfac
2405 ad_ri(i)=0.0_r8
2406
2407# ifdef COOL_SKIN
2408!^ tl_delTc(i)=0.0_r8
2409 ad_deltc(i)=0.0_r8
2410# endif
2411
2412!^ tl_delTc(i)=0.0_r8
2413 ad_deltc(i)=0.0_r8
2414
2415!^ tl_CC(i)=vonKar*tl_Ct(i)/Cd-tl_Cd*CC(i)/Cd
2416 adfac=ad_cc(i)/cd
2417 ad_ct(i)=ad_ct(i)+vonkar*adfac
2418 ad_cd=ad_cd-cc(i)*adfac
2419 ad_cc(i)=0.0_r8
2420
2421 fac=log(blk_zt(ng)/zot10(i))
2422
2423!^ tl_Ct(i)=-tl_fac*Ct(i)/fac
2424 ad_fac=-ad_ct(i)*ct(i)/fac
2425
2426!^ tl_fac=-tl_ZoT10(i)/ZoT10(i)
2427 ad_zot10(i)=ad_zot10(i)-ad_fac/zot10(i)
2428 ad_fac=0.0_r8
2429
2430!
2431! Adjoint Neutral coefficients.
2432!
2433 fac=log(blk_zw(ng)/zo10(i))
2434
2435!^ tl_Cd=-2.0_r8*tl_fac*Cd/fac
2436 ad_fac=-2.0_r8*ad_cd/fac
2437 ad_cd=0.0_r8
2438
2439!^ tl_fac=-tl_Zo10(i)/Zo10(i)
2440 ad_zo10(i)=ad_zo10(i)-ad_fac/zo10(i)
2441 ad_fac=0.0_r8
2442
2443 fac=vonkar/ct10(i)
2444
2445!^ tl_ZoT10(i)=-tl_fac*ZoT10(i)
2446 ad_fac=-ad_zot10(i)*zot10(i)
2447 ad_fac=0.0_r8
2448
2449!^ tl_fac=-tl_Ct10(i)*fac/Ct10(i)
2450 ad_ct10(i)=ad_ct10(i)-fac*ad_fac/ct10(i)
2451 ad_fac=0.0_r8
2452
2453!^ tl_Ct10(i)=-0.5_r8*tl_Cd10(i)*Ct10(i)/Cd10(i)
2454 ad_cd10(i)=ad_cd10(i)-0.5_r8*ct10(i)*ad_ct10(i)/cd10(i)
2455 ad_ct10(i)=0.0_r8
2456
2457 fac=log(10.0_r8/zo10(i))
2458
2459!^ tl_Cd10(i)=-2.0_r8*tl_fac*Cd10(i)/fac
2460 ad_fac=-2.0_r8*cd10(i)*ad_cd10(i)/fac
2461 ad_cd10(i)=0.0_r8
2462
2463!^ tl_fac=-tl_Zo10(i)/Zo10(i)
2464 ad_zo10(i)=ad_zo10(i)-ad_fac/zo10(i)
2465 ad_fac=0.0_r8
2466
2467!^ tl_Zo10(i)=0.022_r8*tl_Wstar(i)*Wstar(i)/g- &
2468!^ & tl_Wstar(i)*0.11_r8*VisAir(i)/(Wstar(i)*Wstar(i))
2469
2470 ad_wstar(i)=ad_wstar(i)+0.022_r8*wstar(i)*ad_zo10(i)/g- &
2471 & 0.11_r8*visair(i)*ad_zo10(i)/(wstar(i)*wstar(i))
2472 ad_zo10(i)=0.0_r8
2473
2474!^ tl_Wstar(i)=0.035_r8*tl_u10(i)
2475 ad_u10(i)=ad_u10(i)+0.035_r8*ad_wstar(i)
2476 ad_wstar(i)=0.0_r8
2477
2478!^ tl_u10(i)=0.0_r8
2479
2480 ad_u10(i)=0.0_r8
2481
2482!^ tl_ZoW(i)=0.0_r8
2483 ad_zow(i)=0.0_r8
2484
2485
2486!
2487! Adjoint wind is measured relative to sea surface and include
2488! gustiness.
2489!
2490
2491!^ tl_delT(i)=tl_TseaC(i)
2492 ad_tseac(i)=ad_tseac(i)+ad_delt(i)
2493 ad_delt(i)=0.0_r8
2494
2495!^ tl_delQ(i)=tl_Qsea(i)
2496 ad_qsea(i)=ad_qsea(i)+ad_delq(i)
2497 ad_delq(i)=0.0_r8
2498
2499!^ tl_delW(i)=0.0_r8
2500 ad_delw(i)=0.0_r8
2501
2502!^ tl_Wgus(i)=0.0_r8
2503 ad_wgus(i)=0.0_r8
2504
2505!
2506! Adjoint latent heat of vaporization (J/kg) at sea surface, Hlv.
2507!
2508!^ tl_Hlv(i,j)=-0.00237_r8*tl_TseaC(i)*1.0E+6_r8
2509!^
2510 ad_tseac(i)=ad_tseac(i)-0.00237_r8*1.0e+6_r8*ad_hlv(i,j)
2511 ad_hlv(i,j)=0.0_r8
2512!
2513! Adjoint Qsea (kg/kg) from vapor pressure.
2514!
2515 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
2516 & exp(17.502_r8*tseac(i)/(240.97_r8+tseac(i)))
2517 cff=cff*0.98_r8
2518!^ tl_Qsea(i)=tl_fac*(1.0_r8+0.378_r8*cff/((PairM-0.378_r8*cff)))
2519!^
2520 ad_fac=ad_qsea(i)*(1.0_r8+0.378_r8*cff/((pairm-0.378_r8*cff)))
2521 ad_qsea(i)=0.0_r8
2522!^ tl_fac=0.62197_r8*(tl_cff/(PairM-0.378_r8*cff))
2523!^
2524 ad_cff=ad_cff+0.62197_r8*ad_fac/(pairm-0.378_r8*cff)
2525 ad_fac=0.0_r8
2526!
2527! Adjoint Vapor Pressure reduced for salinity
2528! (Kraus and Businger, 1994, pp 42).
2529!
2530!^ tl_cff=tl_cff*0.98_r8
2531!^
2532 ad_cff=ad_cff*0.98_r8
2533!
2534! Adjoint water saturation vapor pressure (mb), using Teten formula.
2535!
2536 fac=17.502_r8*tseac(i)/(240.97_r8+tseac(i))
2537 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
2538 & exp(17.502_r8*tseac(i)/(240.97_r8+tseac(i)))
2539!^ tl_cff=tl_fac*cff
2540!^
2541 ad_fac=ad_cff*cff
2542 ad_cff=0.0_r8
2543!^ tl_fac=17.502_r8*tl_TseaC(i)/(240.97_r8+TseaC(i))- &
2544!^ & tl_TseaC(i)*fac/(240.97_r8+TseaC(i))
2545!^
2546 ad_tseac(i)=ad_tseac(i)+ &
2547 & 17.502_r8*ad_fac/(240.97_r8+tseac(i))- &
2548 & fac*ad_fac/(240.97_r8+tseac(i))
2549 ad_fac=0.0_r8
2550!
2551!-----------------------------------------------------------------------
2552! Adjoint net longwave radiation (W/m2), LRad.
2553!-----------------------------------------------------------------------
2554!
2555# ifdef MASKING
2556!^ tl_LRad(i,j)=tl_LRad(i,j)*rmask(i,j)
2557!^
2558 ad_lrad(i,j)=ad_lrad(i,j)*rmask(i,j)
2559# endif
2560# if defined LONGWAVE
2561!
2562! Use Berliand (1952) formula to calculate net longwave radiation.
2563! The equation for saturation vapor pressure is from Gill (Atmosphere-
2564! Ocean Dynamics, pp 606). Here the coefficient in the cloud term
2565! is assumed constant, but it is a function of latitude varying from
2566! 1.0 at poles to 0.5 at the Equator).
2567!
2568 cff2=tairk(i)*tairk(i)*tairk(i)
2569!^ tl_LRad(i,j)=-emmiss*StefBo*cff2*4.0_r8*tl_TseaK(i)
2570!^
2571 ad_tseak(i)=ad_tseak(i)-emmiss*stefbo*cff2*4.0_r8*ad_lrad(i,j)
2572 ad_lrad(i,j)=0.0_r8
2573
2574# elif defined LONGWAVE_OUT
2575!
2576! Treat input longwave data as downwelling radiation only and add
2577! outgoing IR from model sea surface temperature.
2578!
2579!^ tl_LRad(i,j)=tl_lrflx(i,j)*Hscale- &
2580!^ & 4.0_r8*emmiss*StefBo* &
2581!^ & tl_TseaK(i)*TseaK(i)*TseaK(i)*TseaK(i)
2582!^
2583 ad_lrflx(i,j)=ad_lrflx(i,j)+ad_lrad(i,j)*hscale
2584 ad_tseak(i)=ad_tseak(i)- &
2585 & 4.0_r8*emmiss*stefbo*ad_lrad(i,j)* &
2586 & tseak(i)*tseak(i)*tseak(i)
2587 ad_lrad(i,j)=0.0_r8
2588# else
2589!^ tl_LRad(i,j)=tl_lrflx(i,j)*Hscale
2590!^
2591 ad_lrflx(i,j)=ad_lrflx(i,j)+ad_lrad(i,j)*hscale
2592 ad_lrad(i,j)=0.0_r8
2593# endif
2594!
2595! Initialize.
2596!
2597!^ tl_delTc(i)=0.0_r8
2598!^
2599 ad_delqc(i)=0.0_r8
2600
2601!^ tl_delQc(i)=0.0_r8
2602!^
2603 ad_deltc(i)=0.0_r8
2604!^ tl_LHeat(i,j)=tl_lhflx(i,j)*Hscale
2605!^
2606 ad_lhflx(i,j)=ad_lhflx(i,j)+ad_lheat(i,j)*hscale
2607 ad_lheat(i,j)=0.0_r8
2608!^ tl_SHeat(i,j)=tl_shflx(i,j)*Hscale
2609!^
2610 ad_shflx(i,j)=ad_shflx(i,j)+ad_sheat(i,j)*hscale
2611 ad_sheat(i,j)=0.0_r8
2612!^ tl_Taux(i,j)=0.0_r8
2613!^
2614 ad_taux(i,j)=0.0_r8
2615!^ tl_Tauy(i,j)=0.0_r8
2616!^
2617 ad_tauy(i,j)=0.0_r8
2618!^ tl_Scff(i)=tl_beta(i,j)
2619!^
2620 ad_beta(i,j)=ad_beta(i,j)+ad_scff(i)
2621 ad_scff(i)=0.0_r8
2622!^ tl_Tcff(i)=tl_alpha(i,j)
2623!^
2624 ad_alpha(i,j)=ad_alpha(i,j)+ad_tcff(i)
2625 ad_tcff(i)=0.0_r8
2626!^ tl_rhoSea(i)=tl_rho(i,j,N(ng))
2627!^
2628 ad_rho(i,j,n(ng))=ad_rho(i,j,n(ng))+ad_rhosea(i)
2629 ad_rhosea(i)=0.0_r8
2630!^ tl_TseaK(i)=tl_TseaC(i)
2631!^
2632 ad_tseac(i)=ad_tseac(i)+ad_tseak(i)
2633 ad_tseak(i)=0.0_r8
2634!^ tl_TseaC(i)=tl_t(i,j,N(ng),nrhs,itemp)
2635!^
2636 ad_t(i,j,n(ng),nrhs,itemp)=ad_t(i,j,n(ng),nrhs,itemp)+ &
2637 & ad_tseac(i)
2638 ad_tseac(i)=0.0_r8
2639 END DO
2640 END DO
2641!
2642 RETURN
2643 END SUBROUTINE ad_bulk_flux_tile
2644!
2645 FUNCTION ad_bulk_psiu (tl_ZoL, ZoL, pi)
2646!
2647!=======================================================================
2648! !
2649! This function evaluates the stability function for wind speed !
2650! by matching Kansas and free convection forms. The convective !
2651! form follows Fairall et al. (1996) with profile constants from !
2652! Grachev et al. (2000) BLM. The stable form is from Beljaars !
2653! and Holtslag (1991). !
2654! !
2655!=======================================================================
2656!
2657 USE mod_kinds
2658!
2659! Function result
2660!
2661 real(r8) :: ad_bulk_psiu
2662!
2663! Imported variable declarations.
2664!
2665 real(r8), intent(in) :: tl_zol, zol
2666 real(dp), intent(in) :: pi
2667!
2668! Local variable declarations.
2669!
2670 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
2671
2672 real(r8) :: fw, cff, psic, psik, x, y, fac
2673 real(r8) :: tl_fw, tl_cff, tl_psic, tl_psik, tl_x, tl_y, tl_fac
2674!
2675!-----------------------------------------------------------------------
2676! Compute stability function, PSI.
2677!-----------------------------------------------------------------------
2678!
2679! Unstable conditions.
2680!
2681 IF (zol.lt.0.0_r8) THEN
2682 x=(1.0_r8-15.0_r8*zol)**0.25_r8
2683 tl_x=-0.25_r8*15.0_r8*tl_zol/((1.0_r8-15.0_r8*zol)**0.75_r8)
2684 psik=2.0_r8*log(0.5_r8*(1.0_r8+x))+ &
2685 & log(0.5_r8*(1.0_r8+x*x))- &
2686 & 2.0_r8*atan(x)+0.5_r8*pi
2687 tl_psik=tl_x/(0.5_r8*(1.0_r8+x))+ &
2688 & tl_x*x/(0.5_r8*(1.0_r8+x*x))- &
2689 & 2.0_r8*tl_x/(1.0_r8+x*x)
2690!! & 2.0_r8*tl_x/SQRT(1.0_r8+x*x)
2691!
2692! For very unstable conditions, use free-convection (Fairall).
2693!
2694 cff=sqrt(3.0_r8)
2695 y=(1.0_r8-10.15_r8*zol)**r3
2696 tl_y=-r3*10.15_r8*tl_zol*(1.0_r8-10.15_r8*zol)**(r3-1.0_r8)
2697 fac=(1.0_r8+2.0_r8*y)/cff
2698 tl_fac=2.0_r8*tl_y/cff
2699 psic=1.5_r8*log(r3*(1.0_r8+y+y*y))- &
2700 & cff*atan(fac)+pi/cff
2701 tl_psic=tl_y*(1.0_r8+2.0_r8*y)*1.5_r8/(1.0_r8+y+y*y)- &
2702 & cff*tl_fac/(1.0_r8+fac*fac)
2703!! & cff*tl_fac/SQRT(1.0_r8+fac*fac)
2704!
2705! Match Kansas and free-convection forms with weighting Fw.
2706!
2707 cff=zol*zol
2708 tl_cff=2.0_r8*tl_zol*zol
2709 fw=cff/(1.0_r8+cff)
2710 tl_fw=tl_cff/(1.0_r8+cff)-tl_cff*fw*fw/cff
2711 ad_bulk_psiu=(1.0_r8-fw)*tl_psik-tl_fw*psik &
2712 & +tl_fw*psic+fw*tl_psic
2713!
2714! Stable conditions.
2715!
2716 ELSE
2717 cff=min(50.0_r8,0.35_r8*zol)
2718 tl_cff=(0.5_r8-sign(0.5_r8,0.35_r8*zol-50.0))*0.35_r8*tl_zol
2719 fac=exp(cff)
2720 tl_fac=tl_cff*fac
2721 ad_bulk_psiu=-(tl_zol+0.6667_r8*tl_zol/fac- &
2722 & tl_fac*0.6667_r8*(zol-14.28_r8)/(fac*fac))
2723 END IF
2724!
2725 RETURN
2726 END FUNCTION ad_bulk_psiu
2727!
2728 FUNCTION ad_bulk_psit (tl_ZoL, ZoL, pi)
2729!
2730!=======================================================================
2731! !
2732! This function evaluates the stability function for moisture and !
2733! heat by matching Kansas and free convection forms. The convective !
2734! form follows Fairall et al. (1996) with profile constants from !
2735! Grachev et al. (2000) BLM. The stable form is from Beljaars and !
2736! and Holtslag (1991). !
2737!
2738!=======================================================================
2739! !
2740!
2741 USE mod_kinds
2742!
2743! Function result
2744!
2745 real(r8) :: ad_bulk_psit
2746!
2747! Imported variable declarations.
2748!
2749 real(r8), intent(in) :: tl_zol, zol
2750 real(dp), intent(in) :: pi
2751!
2752! Local variable declarations.
2753!
2754 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
2755
2756 real(r8) :: fw, cff, psic, psik, x, y, fac
2757 real(r8) :: tl_fw, tl_cff, tl_psic, tl_psik, tl_x, tl_y, tl_fac
2758!
2759!-----------------------------------------------------------------------
2760! Compute stability function, PSI.
2761!-----------------------------------------------------------------------
2762!
2763! Unstable conditions.
2764!
2765 IF (zol.lt.0.0_r8) THEN
2766 x=(1.0_r8-15.0_r8*zol)**0.5_r8
2767 IF (x.ne.0.0) THEN
2768 tl_x=-0.5_r8*15.0_r8*tl_zol/x
2769 ELSE
2770 tl_x=0.0_r8
2771 END IF
2772 psik=2.0_r8*log(0.5_r8*(1.0_r8+x))
2773 tl_psik=tl_x/(0.5_r8*(1.0_r8+x))
2774!
2775! For very unstable conditions, use free-convection (Fairall).
2776!
2777 cff=sqrt(3.0_r8)
2778 y=(1.0_r8-34.15_r8*zol)**r3
2779 tl_y=-r3*34.15_r8*tl_zol*(1.0_r8-34.15_r8*zol)**(r3-1.0_r8)
2780 fac=(1.0_r8+2.0_r8*y)/cff
2781 tl_fac=2.0_r8*tl_y/cff
2782 psic=1.5_r8*log(r3*(1.0_r8+y+y*y))- &
2783 & cff*atan(fac)+pi/cff
2784 tl_psic=tl_y*(1.0_r8+2.0_r8*y)*1.5_r8/(1.0_r8+y+y*y)- &
2785 & cff*tl_fac/(1.0_r8+fac*fac)
2786!! & cff*tl_fac/SQRT(1.0_r8+fac*fac)
2787!
2788! Match Kansas and free-convection forms with weighting Fw.
2789!
2790 cff=zol*zol
2791 tl_cff=2.0_r8*tl_zol*zol
2792 fw=cff/(1.0_r8+cff)
2793 tl_fw=tl_cff/(1.0_r8+cff)-tl_cff*fw*fw/cff
2794 ad_bulk_psit=(1.0_r8-fw)*tl_psik-tl_fw*psik &
2795 & +tl_fw*psic+fw*tl_psic
2796!
2797! Stable conditions.
2798!
2799 ELSE
2800 cff=min(50.0_r8,0.35_r8*zol)
2801 tl_cff=(0.5_r8-sign(0.5_r8,0.35_r8*zol-50.0_r8))*0.35_r8*tl_zol
2802 fac=exp(cff)
2803 tl_fac=tl_cff*fac
2804 ad_bulk_psit=-(3.0_r8*tl_zol*(1.0_r8+2.0_r8*zol)**0.5_r8+ &
2805 & 0.6667_r8*tl_zol/fac- &
2806 & tl_fac*0.6667_r8*(zol-14.28_r8)/(fac*fac))
2807 END IF
2808!
2809 RETURN
2810 END FUNCTION ad_bulk_psit
2811#endif
2812 END MODULE ad_bulk_flux_mod
ad_bulk_flux_mod
Definition ad_bulk_flux.F:2
ad_bulk_flux_mod::ad_bulk_flux_tile
subroutine ad_bulk_flux_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs nrhs, rmask, umask, vmask, alpha, ad_alpha, beta, ad_beta, rho, ad_rho, t, ad_t, hair, pair, tair, uwind, vwind, awave, pwave, rain, lhflx, ad_lhflx, lrflx, ad_lrflx, shflx, ad_shflx, srflx, stflx, ad_stflx, evap, ad_evap, ad_sustr, ad_svstr)
Definition ad_bulk_flux.F:160
ad_bulk_flux_mod::ad_bulk_psit
real(r8) function ad_bulk_psit(tl_zol, zol, pi)
Definition ad_bulk_flux.F:2729
ad_bulk_flux_mod::ad_bulk_psiu
real(r8) function ad_bulk_psiu(tl_zol, zol, pi)
Definition ad_bulk_flux.F:2646
ad_bulk_flux_mod::ad_bulk_flux
subroutine, public ad_bulk_flux(ng, tile)
Definition ad_bulk_flux.F:55
bulk_flux_mod
Definition bulk_flux.F:2
bulk_flux_mod::bulk_psiu
real(r8) function, public bulk_psiu(zol, pi)
Definition bulk_flux.F:1466
bulk_flux_mod::bulk_psit
real(r8) function, public bulk_psit(zol, pi)
Definition bulk_flux.F:1531
exchange_2d_mod
Definition exchange_2d.F:2
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_kinds
Definition mod_kinds.F:2
mod_kinds::r8
integer, parameter r8
Definition mod_kinds.F:28
mod_kinds::dp
integer, parameter dp
Definition mod_kinds.F:25
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::iadm
integer, parameter iadm
Definition mod_param.F:665
mod_scalars
Definition mod_scalars.F:2
mod_scalars::blk_dter
real(dp) blk_dter
Definition mod_scalars.F:1501
mod_scalars::blk_zt
real(r8), dimension(:), allocatable blk_zt
Definition mod_scalars.F:1506
mod_scalars::blk_zw
real(r8), dimension(:), allocatable blk_zw
Definition mod_scalars.F:1507
mod_scalars::vonkar
real(dp) vonkar
Definition mod_scalars.F:469
mod_scalars::blk_cpw
real(dp) blk_cpw
Definition mod_scalars.F:1497
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::cp
real(dp) cp
Definition mod_scalars.F:456
mod_scalars::blk_zabl
real(dp) blk_zabl
Definition mod_scalars.F:1499
mod_scalars::blk_visw
real(dp) blk_visw
Definition mod_scalars.F:1503
mod_scalars::blk_tcw
real(dp) blk_tcw
Definition mod_scalars.F:1502
mod_scalars::blk_rgas
real(dp) blk_rgas
Definition mod_scalars.F:1498
mod_scalars::isalt
integer isalt
Definition mod_scalars.F:127
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
mod_scalars::stefbo
real(dp) stefbo
Definition mod_scalars.F:460
mod_scalars::rhow
real(dp) rhow
Definition mod_scalars.F:462
mod_scalars::blk_zq
real(r8), dimension(:), allocatable blk_zq
Definition mod_scalars.F:1505
mod_scalars::blk_beta
real(dp) blk_beta
Definition mod_scalars.F:1500
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_scalars::blk_cpa
real(dp) blk_cpa
Definition mod_scalars.F:1496
mod_scalars::emmiss
real(dp) emmiss
Definition mod_scalars.F:461
mod_scalars::pi
real(dp), parameter pi
Definition mod_scalars.F:813
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::ad_mp_exchange2d
subroutine ad_mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, ad_a, ad_b, ad_c, ad_d)
Definition mp_exchange.F:4058
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3