forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
exf_getforcing.F
309 lines (280 loc) · 9.91 KB
/
exf_getforcing.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_getforcing.F,v 1.49 2017/01/27 17:22:55 jmc Exp $
C $Name: $
#include "EXF_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
CBOI
C
C !TITLE: EXTERNAL FORCING
C !AUTHORS: mitgcm developers ( mitgcm-support@mitgcm.org )
C !AFFILIATION: Massachussetts Institute of Technology
C !DATE:
C !INTRODUCTION: External forcing package
C \bv
C * The external forcing package, in conjunction with the
C calendar package (cal), enables the handling of realistic forcing
C fields of differing temporal forcing patterns.
C * It comprises climatological restoring and relaxation
C * Bulk formulae are implemented to convert atmospheric fields
C to surface fluxes.
C * An interpolation routine provides on-the-fly interpolation of
C forcing fields an arbitrary grid onto the model grid.
C * A list of EXF variables and units is in EXF_FIELDS.h
C
C !CALLING SEQUENCE:
C ...
C EXF_GETFORCING (TOP LEVEL ROUTINE)
C |
C |-- EXF_GETCLIM (get climatological fields used e.g. for relax.)
C | |--- exf_set_climtemp (relax. to 3-D temperature field)
C | |--- exf_set_climsalt (relax. to 3-D salinity field)
C | |--- exf_set_climsst (relax. to 2-D SST field)
C | |--- exf_set_climsss (relax. to 2-D SSS field)
C | o
C |
C |-- EXF_GETFFIELDS <- this one does almost everything
C | | 1. reads in fields, either flux or atmos. state,
C | | depending on CPP options (for each variable two fields
C | | consecutive in time are read in and interpolated onto
C | | current time step).
C | | 2. If forcing is atmos. state and control is atmos. state,
C | | then the control variable anomalies are read here
C | | * ctrl_getatemp
C | | * ctrl_getaqh
C | | * ctrl_getuwind
C | | * ctrl_getvwind
C | | If forcing and control are fluxes, then
C | | controls are added later.
C | o
C |
C |-- EXF_CHECK_RANGE
C | | Check whether read fields are within assumed range
C | | (may capture mismatches in units)
C | o
C |
C |-- EXF_RADIATION
C | | Compute net or downwelling radiative fluxes via
C | | Stefan-Boltzmann law in case only one is known.
C |-- EXF_WIND
C | | Compute air-sea wind-stress from winds (or the other way)
C |-- EXF_BULKFORMULAE
C | | Compute air-sea buoyancy fluxes from atmospheric
C | | state following Large and Pond, JPO, 1981/82
C | o
C |
C |-- < add time-mean river runoff here, if available >
C |
C |-- < update tile edges here >
C |
C |-- EXF_GETSURFACEFLUXES
C | | If forcing and control are fluxes, then
C | | control vector anomalies are added here.
C | |--- ctrl_get_gen
C | o
C |
C |-- < treatment of hflux w.r.t. swflux >
C |
C |-- EXF_DIAGNOSTICS_FILL
C | | Do EXF-related diagnostics output here.
C |-- EXF_MONITOR
C | | Monitor EXF-forcing fields
C | o
C |
C |-- EXF_MAPFIELDS
C | | Forcing fields from exf package are mapped onto
C | | mitgcm forcing arrays.
C | | Mapping enables a runtime rescaling of fields
C | o
C
C \ev
CEOI
CBOP
C !ROUTINE: EXF_GETFORCING
C !INTERFACE:
SUBROUTINE EXF_GETFORCING( myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *=================================================================
C | SUBROUTINE EXF_GETFORCING
C *=================================================================
C o Get the forcing fields for the current time step. The switches
C for the inclusion of the individual forcing components have to
C be set in EXF_OPTIONS.h (or ECCO_CPPOPTIONS.h).
C A note on surface fluxes:
C The MITgcm-UV vertical coordinate z is positive upward.
C This implies that a positive flux is out of the ocean
C model. However, the wind stress forcing is not treated
C this way. A positive zonal wind stress accelerates the
C model ocean towards the east.
C started: eckert@mit.edu, heimbach@mit.edu, ralf@ocean.mit.edu
C mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
C *=================================================================
C | SUBROUTINE EXF_GETFORCING
C *=================================================================
C \ev
C !USES:
IMPLICIT NONE
C == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#include "EXF_PARAM.h"
#include "EXF_FIELDS.h"
#include "EXF_CONSTANTS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C == routine arguments ==
_RL myTime
INTEGER myIter
INTEGER myThid
C !LOCAL VARIABLES:
C == local variables ==
INTEGER bi,bj
INTEGER i,j,k
C == end of interface ==
CEOP
C Get values of climatological fields.
CALL EXF_GETCLIM( myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef ALLOW_ATM_TEMP
CADJ STORE precip0 = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE precip1 = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE snowprecip0 = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE snowprecip1 = comlev1, key=ikey_dynamics, kind=isbyte
# endif
#endif
C Get the surface forcing fields.
CALL EXF_GETFFIELDS( myTime, myIter, myThid )
IF ( .NOT.useAtmWind ) THEN
IF ( stressIsOnCgrid .AND. ustressfile.NE.' '
& .AND. vstressfile.NE.' ' )
& CALL EXCH_UV_XY_RL( ustress, vstress, .TRUE., myThid )
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef ALLOW_AUTODIFF_MONITOR
CALL EXF_ADJOINT_SNAPSHOTS( 2, myTime, myIter, myThid )
# endif
#endif
#ifdef ALLOW_DOWNWARD_RADIATION
C Set radiative fluxes
CALL EXF_RADIATION( myTime, myIter, myThid )
#endif
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ustress = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE vstress = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE uwind = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE vwind = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE wspeed = comlev1, key=ikey_dynamics, kind=isbyte
#endif
C Set wind fields
CALL EXF_WIND( myTime, myIter, myThid )
#ifdef ALLOW_ATM_TEMP
# ifdef ALLOW_BULKFORMULAE
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ustress = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE vstress = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE uwind = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE vwind = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE wspeed = comlev1, key=ikey_dynamics, kind=isbyte
# endif
C Compute turbulent fluxes (and surface stress) from bulk formulae
CALL EXF_BULKFORMULAE( myTime, myIter, myThid )
# endif /* ALLOW_BULKFORMULAE */
#endif /* ALLOW_ATM_TEMP */
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
#ifdef ALLOW_ATM_TEMP
C compute hflux & sflux from multiple components
DO j = 1,sNy
DO i = 1,sNx
C Net surface heat flux.
hflux(i,j,bi,bj) =
& - hs(i,j,bi,bj)
& - hl(i,j,bi,bj)
& + lwflux(i,j,bi,bj)
#ifndef SHORTWAVE_HEATING
& + swflux(i,j,bi,bj)
#endif
C fresh-water flux from Precipitation and Evaporation.
sflux(i,j,bi,bj) = evap(i,j,bi,bj) - precip(i,j,bi,bj)
ENDDO
ENDDO
#endif /* ALLOW_ATM_TEMP */
C Apply runoff, masks and exchanges
k = 1
DO j = 1,sNy
DO i = 1,sNx
#ifdef ALLOW_RUNOFF
sflux(i,j,bi,bj) = sflux(i,j,bi,bj) - runoff(i,j,bi,bj)
#endif
hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskC(i,j,k,bi,bj)
sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskC(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
C Update the tile edges: needed for some EXF fields involved in horizontal
C averaging, e.g., wind-stress; fields used by main model or other pkgs
C are exchanged in EXF_MAPFIELDS.
c _EXCH_XY_RL(hflux, myThid)
c _EXCH_XY_RL(sflux, myThid)
IF ( stressIsOnCgrid ) THEN
CALL EXCH_UV_XY_RL( ustress, vstress, .TRUE., myThid )
ELSE
CALL EXCH_UV_AGRID_3D_RL(ustress, vstress, .TRUE., 1, myThid)
ENDIF
#ifdef SHORTWAVE_HEATING
c _EXCH_XY_RL(swflux, myThid)
#endif
#ifdef ATMOSPHERIC_LOADING
c _EXCH_XY_RL(apressure, myThid)
#endif
#ifdef EXF_SEAICE_FRACTION
c _EXCH_XY_RL(areamask, myThid)
#endif
C Get values of the surface flux anomalies.
CALL EXF_GETSURFACEFLUXES( myTime, myIter, myThid )
IF ( useExfCheckRange .AND.
& ( myIter.EQ.nIter0 .OR. exf_debugLev.GE.debLevC ) ) THEN
CALL EXF_CHECK_RANGE( myTime, myIter, myThid )
ENDIF
#ifdef ALLOW_AUTODIFF
# ifdef ALLOW_AUTODIFF_MONITOR
CALL EXF_ADJOINT_SNAPSHOTS( 1, myTime, myIter, myThid )
# endif
#endif /* ALLOW_AUTODIFF */
#ifdef ALLOW_ATM_TEMP
# ifdef SHORTWAVE_HEATING
C Treatment of qnet
C The location of the summation of Qnet in exf_mapfields is unfortunate:
C For backward compatibility issues we want it to happen after
C applying control variables, but before exf_diagnostics_fill.
C Therefore, we DO it exactly here:
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j = 1-oLy,sNy+oLy
DO i = 1-oLx,sNx+oLx
hflux(i,j,bi,bj) = hflux(i,j,bi,bj) + swflux(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
# endif /* SHORTWAVE_HEATING */
#endif /* ALLOW_ATM_TEMP */
C Diagnostics output
CALL EXF_DIAGNOSTICS_FILL( myTime, myIter, myThid )
C Monitor output
CALL EXF_MONITOR( myTime, myIter, myThid )
C Map the forcing fields onto the corresponding model fields.
CALL EXF_MAPFIELDS( myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF
# ifdef ALLOW_AUTODIFF_MONITOR
IF ( .NOT. useSEAICE )
& CALL EXF_ADJOINT_SNAPSHOTS( 3, myTime, myIter, myThid )
# endif
#endif /* ALLOW_AUTODIFF */
RETURN
END