1 | PROGRAM main_bpc_hpg |
---|
2 | |
---|
3 | !-------------------------------------------------------------------------------------------- |
---|
4 | ! 1. Declarations |
---|
5 | |
---|
6 | USE par_kind |
---|
7 | USE netcdf |
---|
8 | USE len_oce |
---|
9 | USE phycst |
---|
10 | USE in_out_manager ! I/O manager |
---|
11 | USE eosinsitu |
---|
12 | USE zpshde |
---|
13 | USE dynhpgs |
---|
14 | USE zdfmxl |
---|
15 | USE ldfslp |
---|
16 | USE traldf_iso |
---|
17 | USE traldf_iso_tile |
---|
18 | USE tiletype |
---|
19 | |
---|
20 | IMPLICIT NONE |
---|
21 | |
---|
22 | INTEGER :: ji, jj, jk, jt |
---|
23 | INTEGER :: jpi_in, jpj_in, jpk_in |
---|
24 | INTEGER :: istep,maxstep |
---|
25 | INTEGER :: iostat, ncid, dimid, varid, chunksize |
---|
26 | INTEGER :: x_dimid, y_dimid, z_dimid, i_dimid, h_dimid |
---|
27 | INTEGER :: x_varid, y_varid, z_varid, nav_lon_varid, nav_lat_varid, nav_lev_varid |
---|
28 | INTEGER :: hpgi_varid, hpgj_varid, ta_varid, sa_varid, hbpcgi_varid, hbpcgj_varid |
---|
29 | INTEGER :: numnam = 12 ! unit number for namelist |
---|
30 | INTEGER :: kpass = 1 ! for simple harmonic isopycnal diffusion |
---|
31 | INTEGER :: nlb10 = 7 ! w index just below 10 m depth |
---|
32 | |
---|
33 | REAL :: fill = -99999.0 ! note that using REAL(wp) here gave errors in NF90_ENDDEF |
---|
34 | REAL(wp) :: r1_2 = 0.5_wp |
---|
35 | REAL(wp) :: r2dt = 7200. ! 1 hour timestep in seconds |
---|
36 | REAL(wp) :: zcoef, zUfac |
---|
37 | |
---|
38 | INTEGER :: jsi, jsj, jsk ! tile loop indices (copied into tile) |
---|
39 | TYPE(tile_type) :: tile |
---|
40 | |
---|
41 | ! Fields read in from files |
---|
42 | |
---|
43 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: in3d |
---|
44 | REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: in4d |
---|
45 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: nav_lon, nav_lat |
---|
46 | REAL(wp), DIMENSION(:), ALLOCATABLE :: nav_lev |
---|
47 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e1t, e1u, e1v, e2t, e2u, e2v |
---|
48 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3t_0, e3u_0, e3v_0, e3w_0 |
---|
49 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmask, umask, vmask, wmask |
---|
50 | REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ts_n, ts_pcbias |
---|
51 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3t_n |
---|
52 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sshn |
---|
53 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: hpgi, hpgj, hbpcgi, hbpcgj |
---|
54 | |
---|
55 | ! Fields calculated after input files have been read in |
---|
56 | |
---|
57 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r1_e1u, r1_e2v, e1e2t, r1_e1e2t, e1e2u, e1e2v, r1_e1e2u, r1_e1e2v, e1_e2v, e2_e1u |
---|
58 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3u_n, e3v_n, e3w_n, gdept_n, gdepw_n, gde3w_n, rhd, rhop, rn2b, avt, ahtu, ahtv |
---|
59 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: gtsu, gtsv, gtui, gtvi |
---|
60 | INTEGER, DIMENSION(:,:), ALLOCATABLE :: mbkt, mbku, mbkv, mikt, miku, mikv, nmln |
---|
61 | REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: ssmask, gru, grv, risfdep, hmld, hmlp, hmlpt, uslpml, vslpml, wslpiml, wslpjml |
---|
62 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uslp, vslp, wslpi, wslpj, akz, ah_wslp2, omlmask |
---|
63 | REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ts_b, ts_a, rab_b |
---|
64 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ua, va |
---|
65 | |
---|
66 | REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zftu, zftv ! arrays used for diagnostics |
---|
67 | REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos |
---|
68 | REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zti, ztj ! |
---|
69 | REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zhpi, zhpj |
---|
70 | |
---|
71 | CHARACTER(LEN=256) :: meshfile, modelfile, pcbiasfile, outfile, out_diag |
---|
72 | |
---|
73 | LOGICAL :: ln_linssh ! True => linear free surface. e3t does not change with time |
---|
74 | LOGICAL :: ln_zco, ln_zps, ln_sco ! choice of hpg scheme. Only one logical should be true ! |
---|
75 | REAL(wp):: rn_avt0, rn_Ud |
---|
76 | |
---|
77 | LOGICAL :: ln_isfcav ! set to false |
---|
78 | LOGICAL :: ln_traldf_msc, ln_traldf_blp, ln_traldf_lap |
---|
79 | |
---|
80 | LOGICAL :: l_hst, l_ptr ! logicals for diagnostics |
---|
81 | |
---|
82 | INTEGER :: ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile |
---|
83 | |
---|
84 | LOGICAL :: ln_tiled |
---|
85 | |
---|
86 | REAL(wp) :: sum_ua,sum_va,sum_ts_a, step_time |
---|
87 | REAL(wp) :: et |
---|
88 | |
---|
89 | !-------------------------------------------------------------------------------------------- |
---|
90 | ! 2. User inputs |
---|
91 | ! They include the namelist nam_pcb_hpg |
---|
92 | ! the namelist nameos in eos_init |
---|
93 | ! |
---|
94 | |
---|
95 | NAMELIST/nam_bpc_hpg/ meshfile, modelfile, pcbiasfile, outfile, out_diag, ln_linssh, ln_zco, ln_zps, ln_sco, jpi_in, jpj_in, jpk_in, & |
---|
96 | & rn_avt0, rn_Ud, ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile, & |
---|
97 | & maxstep,ln_tiled |
---|
98 | |
---|
99 | write(0,*) 'Starting' |
---|
100 | |
---|
101 | meshfile = 'mesh_mask.nc' |
---|
102 | modelfile = 'restart.nc' |
---|
103 | pcbiasfile = 'pcbias.nc' |
---|
104 | outfile = 'hpg_Jan_2019.nc' |
---|
105 | out_diag = 'out_diag_Jan_2019.txt' |
---|
106 | ln_linssh = .False. |
---|
107 | ln_zco = .False. ; ln_zps = .False. ; ln_sco = .False. |
---|
108 | |
---|
109 | ji_len_3D_tile = 128 ; jj_len_3D_tile = 64 ; jk_len_3D_tile = 3 ! default values for tile lengths |
---|
110 | |
---|
111 | rn_avt0 = 1.2E-5 ! values used in experiment u-bd189 (need to check what this is !!!) |
---|
112 | rn_ud = 0.011 ! values used in experiment u-bd189 |
---|
113 | |
---|
114 | ln_isfcav = .False. |
---|
115 | |
---|
116 | ln_traldf_msc = .False.; ln_traldf_blp = .False.; ln_traldf_lap = .True. |
---|
117 | |
---|
118 | l_hst = .False. ; l_ptr = .False. |
---|
119 | |
---|
120 | nlb10 = 7 |
---|
121 | |
---|
122 | ln_tiled=.false. |
---|
123 | maxstep=10 |
---|
124 | |
---|
125 | OPEN( UNIT=numnam, FILE='nam_bpc_hpg', FORM='FORMATTED', STATUS='OLD' ) |
---|
126 | READ( numnam, nam_bpc_hpg ) |
---|
127 | CLOSE( numnam ) |
---|
128 | |
---|
129 | OPEN( UNIT=6, FILE=out_diag, FORM='FORMATTED', STATUS='UNKNOWN' ) |
---|
130 | WRITE(6, nam_bpc_hpg ) |
---|
131 | |
---|
132 | !-------------------------------------------------------------------------------------------- |
---|
133 | ! 3. Read in the meshmask data |
---|
134 | |
---|
135 | chunksize = 32000000 |
---|
136 | |
---|
137 | write(0,*) 'Opening mesh' |
---|
138 | |
---|
139 | iostat = NF90_OPEN(TRIM(meshfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize) |
---|
140 | |
---|
141 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening meshfile: ', TRIM(nf90_strerror(iostat)) |
---|
142 | |
---|
143 | iostat = NF90_INQ_DIMID(ncid, 'x', dimid) |
---|
144 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpi) |
---|
145 | iostat = NF90_INQ_DIMID(ncid, 'y', dimid) |
---|
146 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpj) |
---|
147 | iostat = NF90_INQ_DIMID(ncid, 'z', dimid) |
---|
148 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpk) |
---|
149 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading dimensions in meshfile: ', TRIM(nf90_strerror(iostat)) |
---|
150 | |
---|
151 | write(0,*) 'jpi, jpj, jpk = ', jpi, jpj, jpk |
---|
152 | |
---|
153 | IF ( jpi /= jpi_in .OR. jpj /= jpj_in .OR. jpk /= jpk_in ) THEN |
---|
154 | write(0,*) ' FATAL ERROR: the input dimensions in the meshmask file and namelist do not match.' |
---|
155 | write(0,*) ' Change the nam_pcb_hpg namelist or the input file ' |
---|
156 | STOP |
---|
157 | END IF |
---|
158 | |
---|
159 | jpim1 = jpi - 1 |
---|
160 | jpjm1 = jpj - 1 |
---|
161 | jpkm1 = jpk - 1 |
---|
162 | |
---|
163 | ALLOCATE( in3d(jpi,jpj,1), in4d(jpi,jpj,jpk,1) ) |
---|
164 | |
---|
165 | ALLOCATE( nav_lon(jpi,jpj), nav_lat(jpi,jpj), nav_lev(jpk), & |
---|
166 | & e1t(jpi,jpj), e1u(jpi,jpj), e1v(jpi,jpj), e2t(jpi,jpj), e2u(jpi,jpj), e2v(jpi,jpj) ) |
---|
167 | |
---|
168 | ALLOCATE( e3t_0(jpi,jpj,jpk), e3u_0(jpi,jpj,jpk), e3v_0(jpi,jpj,jpk), e3w_0(jpi,jpj,jpk), & |
---|
169 | & tmask(jpi,jpj,jpk), umask(jpi,jpj,jpk), vmask(jpi,jpj,jpk), wmask(jpi,jpj,jpk) ) |
---|
170 | |
---|
171 | ALLOCATE( ts_n(jpi,jpj,jpk,jpts), ts_pcbias(jpi,jpj,jpk,jpts) ) |
---|
172 | ALLOCATE( e3t_n(jpi,jpj,jpk), sshn(jpi,jpj) ) |
---|
173 | |
---|
174 | ALLOCATE( hpgi(jpi,jpj,jpk), hpgj(jpi,jpj,jpk), hbpcgi(jpi,jpj,jpk), hbpcgj(jpi,jpj,jpk) ) |
---|
175 | |
---|
176 | ALLOCATE( zri(jpi,jpj), zrj(jpi,jpj), zhi(jpi,jpj), zhj(jpi,jpj), zti(jpi,jpj,jpts), ztj(jpi,jpj,jpts), & |
---|
177 | & zhpi(jpi,jpj,jpk), zhpj(jpi,jpj,jpk)) |
---|
178 | |
---|
179 | iostat = NF90_INQ_VARID(ncid, 'nav_lon', varid) |
---|
180 | iostat = NF90_GET_VAR(ncid, varid, nav_lon(:,:)) |
---|
181 | |
---|
182 | iostat = NF90_INQ_VARID(ncid, 'nav_lat', varid) |
---|
183 | iostat = NF90_GET_VAR(ncid, varid, nav_lat(:,:)) |
---|
184 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading nav_lon or nav_lat: ', TRIM(nf90_strerror(iostat)) |
---|
185 | |
---|
186 | iostat = NF90_INQ_VARID(ncid, 'nav_lev', varid) |
---|
187 | iostat = NF90_GET_VAR(ncid, varid, nav_lev(:)) |
---|
188 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading nav_lev: ', TRIM(nf90_strerror(iostat)) |
---|
189 | |
---|
190 | iostat = NF90_INQ_VARID(ncid, 'e1t', varid) |
---|
191 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
---|
192 | !$OMP PARALLEL WORKSHARE |
---|
193 | e1t(:,:) = in3d(:,:,1) |
---|
194 | !$OMP END PARALLEL WORKSHARE |
---|
195 | |
---|
196 | iostat = NF90_INQ_VARID(ncid, 'e1u', varid) |
---|
197 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
---|
198 | !$OMP PARALLEL WORKSHARE |
---|
199 | e1u(:,:) = in3d(:,:,1) |
---|
200 | !$OMP END PARALLEL WORKSHARE |
---|
201 | |
---|
202 | iostat = NF90_INQ_VARID(ncid, 'e1v', varid) |
---|
203 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
---|
204 | !$OMP PARALLEL WORKSHARE |
---|
205 | e1v(:,:) = in3d(:,:,1) |
---|
206 | !$OMP END PARALLEL WORKSHARE |
---|
207 | |
---|
208 | iostat = NF90_INQ_VARID(ncid, 'e2t', varid) |
---|
209 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
---|
210 | !$OMP PARALLEL WORKSHARE |
---|
211 | e2t(:,:) = in3d(:,:,1) |
---|
212 | !$OMP END PARALLEL WORKSHARE |
---|
213 | |
---|
214 | iostat = NF90_INQ_VARID(ncid, 'e2u', varid) |
---|
215 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
---|
216 | !$OMP PARALLEL WORKSHARE |
---|
217 | e2u(:,:) = in3d(:,:,1) |
---|
218 | !$OMP END PARALLEL WORKSHARE |
---|
219 | |
---|
220 | iostat = NF90_INQ_VARID(ncid, 'e2v', varid) |
---|
221 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
---|
222 | !$OMP PARALLEL WORKSHARE |
---|
223 | e2v(:,:) = in3d(:,:,1) |
---|
224 | !$OMP END PARALLEL WORKSHARE |
---|
225 | |
---|
226 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e1, e2 fields: ', TRIM(nf90_strerror(iostat)) |
---|
227 | |
---|
228 | iostat = NF90_INQ_VARID(ncid, 'e3t_0', varid) |
---|
229 | iostat = NF90_GET_VAR(ncid, varid, e3t_0(:,:,:)) |
---|
230 | |
---|
231 | iostat = NF90_INQ_VARID(ncid, 'e3u_0', varid) |
---|
232 | iostat = NF90_GET_VAR(ncid, varid, e3u_0(:,:,:)) |
---|
233 | |
---|
234 | iostat = NF90_INQ_VARID(ncid, 'e3v_0', varid) |
---|
235 | iostat = NF90_GET_VAR(ncid, varid, e3v_0(:,:,:)) |
---|
236 | |
---|
237 | iostat = NF90_INQ_VARID(ncid, 'e3w_0', varid) |
---|
238 | iostat = NF90_GET_VAR(ncid, varid, e3w_0(:,:,:)) |
---|
239 | |
---|
240 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e3 fields: ', TRIM(nf90_strerror(iostat)) |
---|
241 | |
---|
242 | iostat = NF90_INQ_VARID(ncid, 'tmask', varid) |
---|
243 | iostat = NF90_GET_VAR(ncid, varid, tmask(:,:,:)) |
---|
244 | |
---|
245 | iostat = NF90_INQ_VARID(ncid, 'umask', varid) |
---|
246 | iostat = NF90_GET_VAR(ncid, varid, umask(:,:,:)) |
---|
247 | |
---|
248 | iostat = NF90_INQ_VARID(ncid, 'vmask', varid) |
---|
249 | iostat = NF90_GET_VAR(ncid, varid, vmask(:,:,:)) |
---|
250 | |
---|
251 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading mask fields: ', TRIM(nf90_strerror(iostat)) |
---|
252 | |
---|
253 | iostat = NF90_CLOSE(ncid) |
---|
254 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing mesh file: ', TRIM(nf90_strerror(iostat)) |
---|
255 | |
---|
256 | !-------------------------------------------------------------------------------------------- |
---|
257 | ! 4. Read in the models fields required |
---|
258 | |
---|
259 | write(0,*) 'Opening model fields' |
---|
260 | |
---|
261 | iostat = NF90_OPEN(TRIM(modelfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize) |
---|
262 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening main fields file: ', TRIM(nf90_strerror(iostat)) |
---|
263 | |
---|
264 | iostat = NF90_INQ_DIMID(ncid, 'x', dimid) |
---|
265 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpi) |
---|
266 | iostat = NF90_INQ_DIMID(ncid, 'y', dimid) |
---|
267 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpj) |
---|
268 | iostat = NF90_INQ_DIMID(ncid, 'z', dimid) |
---|
269 | iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpk) |
---|
270 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading main field dimensions: ', TRIM(nf90_strerror(iostat)) |
---|
271 | |
---|
272 | write(0,*) 'jpi, jpj, jpk = ', jpi, jpj, jpk |
---|
273 | |
---|
274 | IF ( jpi /= jpi_in .OR. jpj /= jpj_in .OR. jpk /= jpk_in ) THEN |
---|
275 | write(0,*) ' FATAL ERROR: the input dimensions in the model fields file and namelist do not match.' |
---|
276 | write(0,*) ' Change the nam_pcb_hpg namelist or the input file ' |
---|
277 | STOP |
---|
278 | END IF |
---|
279 | |
---|
280 | iostat = NF90_INQ_VARID(ncid, 'tn', varid) |
---|
281 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
---|
282 | !$OMP PARALLEL WORKSHARE |
---|
283 | ts_n(:,:,:,1) = in4d(:,:,:,1) |
---|
284 | !$OMP END PARALLEL WORKSHARE |
---|
285 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_n(:,:,:,1) field: ', TRIM(nf90_strerror(iostat)) |
---|
286 | |
---|
287 | iostat = NF90_INQ_VARID(ncid, 'sn', varid) |
---|
288 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
---|
289 | !$OMP PARALLEL WORKSHARE |
---|
290 | ts_n(:,:,:,2) = in4d(:,:,:,1) |
---|
291 | !$OMP END PARALLEL WORKSHARE |
---|
292 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_n(:,:,:,2) field: ', TRIM(nf90_strerror(iostat)) |
---|
293 | |
---|
294 | IF ( .NOT.(ln_linssh) ) THEN |
---|
295 | iostat = NF90_INQ_VARID(ncid, 'e3t_n', varid) |
---|
296 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
---|
297 | !$OMP PARALLEL WORKSHARE |
---|
298 | e3t_n(:,:,:) = in4d(:,:,:,1) |
---|
299 | !$OMP END PARALLEL WORKSHARE |
---|
300 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e3t_n field: ', TRIM(nf90_strerror(iostat)) |
---|
301 | ELSE |
---|
302 | !$OMP PARALLEL WORKSHARE |
---|
303 | e3t_n(:,:,:) = e3t_0(:,:,:) |
---|
304 | !$OMP END PARALLEL WORKSHARE |
---|
305 | END IF |
---|
306 | |
---|
307 | IF ( ln_sco ) THEN |
---|
308 | iostat = NF90_INQ_VARID(ncid, 'sshn', varid) |
---|
309 | iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:)) |
---|
310 | !$OMP PARALLEL WORKSHARE |
---|
311 | sshn(:,:) = in3d(:,:,1) |
---|
312 | !$OMP END PARALLEL WORKSHARE |
---|
313 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading sshn field: ', TRIM(nf90_strerror(iostat)) |
---|
314 | ELSE |
---|
315 | !$OMP PARALLEL WORKSHARE |
---|
316 | sshn(:,:) = 0._wp |
---|
317 | !$OMP END PARALLEL WORKSHARE |
---|
318 | END IF |
---|
319 | |
---|
320 | iostat = NF90_CLOSE(ncid) |
---|
321 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing main fields (file: ', TRIM(nf90_strerror(iostat)) |
---|
322 | |
---|
323 | write(0,*) 'Opening pcbias fields' |
---|
324 | |
---|
325 | iostat = NF90_OPEN(TRIM(pcbiasfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize) |
---|
326 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening pcbias (file: ', TRIM(nf90_strerror(iostat)) |
---|
327 | |
---|
328 | iostat = NF90_INQ_VARID(ncid, 'tbias_p', varid) |
---|
329 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
---|
330 | !$OMP PARALLEL WORKSHARE |
---|
331 | ts_pcbias(:,:,:,1) = in4d(:,:,:,1) |
---|
332 | !$OMP END PARALLEL WORKSHARE |
---|
333 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_pcbias(:,:,:,1) field: ', TRIM(nf90_strerror(iostat)) |
---|
334 | |
---|
335 | iostat = NF90_INQ_VARID(ncid, 'sbias_p', varid) |
---|
336 | iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:)) |
---|
337 | !$OMP PARALLEL WORKSHARE |
---|
338 | ts_pcbias(:,:,:,2) = in4d(:,:,:,1) |
---|
339 | !$OMP END PARALLEL WORKSHARE |
---|
340 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_pcbias(:,:,:,2) field: ', TRIM(nf90_strerror(iostat)) |
---|
341 | |
---|
342 | iostat = NF90_CLOSE(ncid) |
---|
343 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing pcbias (file: ', TRIM(nf90_strerror(iostat)) |
---|
344 | |
---|
345 | !-------------------------------------------------------------------------------------------- |
---|
346 | ! 5. Open the output file and define the dimensions and output names |
---|
347 | iostat = NF90_CREATE(TRIM(outfile), NF90_CLOBBER, ncid) |
---|
348 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR creating file: ', TRIM(nf90_strerror(iostat)) |
---|
349 | |
---|
350 | iostat = NF90_DEF_DIM(ncid, 'x', jpi, x_dimid) |
---|
351 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat)) |
---|
352 | iostat = NF90_DEF_DIM(ncid, 'y', jpj, y_dimid) |
---|
353 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat)) |
---|
354 | iostat = NF90_DEF_DIM(ncid, 'z', jpk, z_dimid) |
---|
355 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat)) |
---|
356 | |
---|
357 | iostat = NF90_DEF_VAR(ncid, 'nav_lon', NF90_FLOAT, (/ x_dimid,y_dimid /), nav_lon_varid) |
---|
358 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat)) |
---|
359 | iostat = NF90_DEF_VAR(ncid, 'nav_lat', NF90_FLOAT, (/ x_dimid,y_dimid /), nav_lat_varid) |
---|
360 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat)) |
---|
361 | iostat = NF90_DEF_VAR(ncid, 'nav_lev', NF90_FLOAT, (/ z_dimid /), nav_lev_varid) |
---|
362 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat)) |
---|
363 | |
---|
364 | iostat = NF90_DEF_VAR(ncid, 'hpgi', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hpgi_varid) |
---|
365 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hpgi: ', TRIM(nf90_strerror(iostat)) |
---|
366 | iostat = NF90_DEF_VAR(ncid, 'hpgj', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hpgj_varid) |
---|
367 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hpgj: ', TRIM(nf90_strerror(iostat)) |
---|
368 | iostat = NF90_DEF_VAR(ncid, 'hbpcgi', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hbpcgi_varid) |
---|
369 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hbpcgi: ', TRIM(nf90_strerror(iostat)) |
---|
370 | iostat = NF90_DEF_VAR(ncid, 'hbpcgj', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hbpcgj_varid) |
---|
371 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hbpcgj: ', TRIM(nf90_strerror(iostat)) |
---|
372 | |
---|
373 | iostat = NF90_PUT_ATT(ncid, hpgi_varid, '_FillValue', fill) |
---|
374 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining fill hpgi_varid : ', TRIM(nf90_strerror(iostat)) |
---|
375 | iostat = NF90_PUT_ATT(ncid, hpgj_varid, '_FillValue', fill) |
---|
376 | iostat = NF90_PUT_ATT(ncid, hbpcgi_varid, '_FillValue', fill) |
---|
377 | iostat = NF90_PUT_ATT(ncid, hbpcgj_varid, '_FillValue', fill) |
---|
378 | IF(ln_tiled) THEN |
---|
379 | !iostat = NF90_DEF_VAR(ncid, 'ta', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), ta_varid) |
---|
380 | !iostat = NF90_DEF_VAR(ncid, 'sa', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), sa_varid) |
---|
381 | !iostat = NF90_PUT_ATT(ncid, ta_varid, '_FillValue', fill) |
---|
382 | !iostat = NF90_PUT_ATT(ncid, sa_varid, '_FillValue', fill) |
---|
383 | ENDIF |
---|
384 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining fill hbpcgj_varid : ', TRIM(nf90_strerror(iostat)) |
---|
385 | |
---|
386 | iostat = NF90_ENDDEF(ncid) |
---|
387 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR leaving define mode: ', TRIM(nf90_strerror(iostat)) |
---|
388 | !set timers |
---|
389 | hpg_zco_time = 0. |
---|
390 | hpg_zps_time = 0. |
---|
391 | hpg_sco_time = 0. |
---|
392 | zps_hde_time = 0. |
---|
393 | eos_time = 0. |
---|
394 | eos2d_time = 0. |
---|
395 | diffusion_time = 0. |
---|
396 | step_time = 0. |
---|
397 | !-------------------------------------------------------------------------------------------- |
---|
398 | ! 6. Additional allocations for variables calculated in later sections |
---|
399 | |
---|
400 | #ifndef MANAGE |
---|
401 | !$ACC ENTER DATA COPYIN(nav_lev, nav_lon, nav_lat, e1t, e1u, e1v, e2t, e2u, e2v,& |
---|
402 | !$ACC r1_e1u, r1_e2v, e3t_0, e3w_0, tmask, umask, vmask, wmask, sshn, ts_n, ts_pcbias) |
---|
403 | |
---|
404 | !$ACC ENTER DATA CREATE(zti, zhi, zri, ztj, zhj, zrj, zhpi, zhpj, & |
---|
405 | !$ACC rhd, rhop, gtsu, gtsv, mbku, mbkv, gru, grv, gdepw_n, gde3w_n,& |
---|
406 | !$ACC ua, va, hpgi, hpgj, hbpcgi, hbpcgj, e3t_n, e3w_n, gdept_n, r1_e1u, r1_e2v) |
---|
407 | #endif |
---|
408 | |
---|
409 | ALLOCATE( r1_e1u(jpi,jpj), r1_e2v(jpi,jpj), e1e2t(jpi,jpj), r1_e1e2t(jpi,jpj), e1e2u(jpi,jpj), & |
---|
410 | & e1e2v(jpi,jpj), r1_e1e2u(jpi,jpj), r1_e1e2v(jpi,jpj), e1_e2v(jpi,jpj), e2_e1u(jpi,jpj) ) |
---|
411 | |
---|
412 | ALLOCATE( e3u_n(jpi,jpj,jpk), e3v_n(jpi,jpj,jpk), e3w_n(jpi,jpj,jpk) ) |
---|
413 | |
---|
414 | ALLOCATE( gdept_n(jpi,jpj,jpk), gdepw_n(jpi,jpj,jpk), gde3w_n(jpi,jpj,jpk) ) |
---|
415 | |
---|
416 | ALLOCATE( rhd(jpi,jpj,jpk), rhop(jpi,jpj,jpk), rn2b(jpi,jpj,jpk), avt(jpi,jpj,jpk), ahtu(jpi,jpj,jpk), ahtv(jpi,jpj,jpk) ) |
---|
417 | |
---|
418 | ALLOCATE( gtsu(jpi,jpj, jpts), gtsv(jpi,jpj, jpts), gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts) ) |
---|
419 | |
---|
420 | ALLOCATE( mbkt(jpi,jpj), mbku(jpi,jpj), mbkv(jpi,jpj), mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), nmln(jpi,jpj) ) |
---|
421 | |
---|
422 | ALLOCATE( ssmask(jpi,jpj), gru(jpi,jpj), grv(jpi,jpj), risfdep(jpi,jpj), & |
---|
423 | hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), & |
---|
424 | & uslpml(jpi,jpj), vslpml(jpi,jpj), wslpiml(jpi,jpj), wslpjml(jpi,jpj) ) |
---|
425 | |
---|
426 | ALLOCATE( uslp(jpi,jpj,jpk), vslp(jpi,jpj,jpk), wslpi(jpi,jpj,jpk), wslpj(jpi,jpj,jpk) ) |
---|
427 | |
---|
428 | ALLOCATE( akz(jpi,jpj,jpk), ah_wslp2(jpi,jpj,jpk), omlmask(jpi,jpj,jpk) ) |
---|
429 | |
---|
430 | ALLOCATE( ts_b(jpi,jpj,jpk,jpts), ts_a(jpi,jpj,jpk,jpts), rab_b(jpi,jpj,jpk,jpts) ) |
---|
431 | |
---|
432 | ALLOCATE( ua(jpi,jpj,jpk), va(jpi,jpj,jpk) ) |
---|
433 | |
---|
434 | IF(ln_tiled) ALLOCATE( zftu(jpi,jpj,jpk,jpts), zftv(jpi,jpj,jpk,jpts) ) ! arrays used for diagnostics ! tiled |
---|
435 | |
---|
436 | !-------------------------------------------------------------------------------------------- |
---|
437 | ! 6. Calculate mbkt, mbku, mbkv, ice-shelf fields, e3t_n, e3w_n, gdept, r1_e1u, r1_e2v |
---|
438 | |
---|
439 | !$ACC KERNELS |
---|
440 | !$OMP PARALLEL |
---|
441 | !$OMP WORKSHARE |
---|
442 | mbkt(:,:) = 0 ; mbku(:,:) = 0 ; mbkv(:,:) = 0 |
---|
443 | !$OMP END WORKSHARE |
---|
444 | !$OMP DO |
---|
445 | DO jk = 1, jpk |
---|
446 | DO jj = 1, jpj |
---|
447 | DO ji = 1, jpi |
---|
448 | IF ( tmask(ji,jj,jk) == 1._wp ) mbkt(ji,jj) = jk |
---|
449 | IF ( umask(ji,jj,jk) == 1._wp ) mbku(ji,jj) = jk |
---|
450 | IF ( vmask(ji,jj,jk) == 1._wp ) mbkv(ji,jj) = jk |
---|
451 | END DO |
---|
452 | END DO |
---|
453 | END DO |
---|
454 | !$OMP END DO |
---|
455 | |
---|
456 | !$OMP DO |
---|
457 | DO jj = 1, jpj |
---|
458 | DO ji = 1, jpi |
---|
459 | mbkt(ji,jj) = MAX( mbkt(ji,jj) , 1) |
---|
460 | mbku(ji,jj) = MAX( mbku(ji,jj) , 1) |
---|
461 | mbkv(ji,jj) = MAX( mbkv(ji,jj) , 1) |
---|
462 | END DO |
---|
463 | END DO |
---|
464 | !$OMP END DO |
---|
465 | |
---|
466 | ! set ice-shelf fields for no ice ! |
---|
467 | mikt(:,:) = 1 ; miku(:,:) = 1 ; mikv(:,:) = 1 ; risfdep(:,:) = 0._wp |
---|
468 | gtui(:,:,:) = 0._wp ; gtvi(:,:,:) = 0._wp |
---|
469 | |
---|
470 | ! The following code is taken from dommsk: dom_msk |
---|
471 | !$OMP WORKSHARE |
---|
472 | wmask (:,:,1) = tmask(:,:,1) |
---|
473 | !$OMP END WORKSHARE |
---|
474 | !$OMP DO |
---|
475 | DO jk = 2, jpk ! interior values |
---|
476 | wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) |
---|
477 | END DO |
---|
478 | !$OMP END DO |
---|
479 | |
---|
480 | ! The following code is taken from domvvl: dom_vvl_interpol for case 'W' interpolation |
---|
481 | |
---|
482 | IF ( ln_linssh ) THEN |
---|
483 | !$OMP WORKSHARE |
---|
484 | e3w_n(:,:,:) = e3w_0(:,:,:) |
---|
485 | !$OMP END WORKSHARE |
---|
486 | |
---|
487 | ELSE |
---|
488 | !$OMP WORKSHARE |
---|
489 | e3w_n(:,:,1) = e3w_0(:,:,1) + e3t_n(:,:,1) - e3t_0(:,:,1) |
---|
490 | !$OMP END WORKSHARE |
---|
491 | ! - ML - The use of mask in this formula enables the special treatment of the last w-point without indirect adressing |
---|
492 | !$OMP DO |
---|
493 | DO jk = 2, jpk |
---|
494 | e3w_n(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( e3t_n(:,:,jk-1) - e3t_0(:,:,jk-1) ) & |
---|
495 | & + 0.5_wp * tmask(:,:,jk) * ( e3t_n(:,:,jk ) - e3t_0(:,:,jk ) ) |
---|
496 | END DO |
---|
497 | !$OMP END DO |
---|
498 | |
---|
499 | END IF |
---|
500 | |
---|
501 | ! The following code is taken from domhgr: SUBROUTINE dom_hgr |
---|
502 | !$OMP WORKSHARE |
---|
503 | e1e2t (:,:) = e1t(:,:) * e2t(:,:) ; r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) |
---|
504 | e1e2u (:,:) = e1u(:,:) * e2u(:,:) ; r1_e1e2u(:,:) = 1._wp / e1e2u(:,:) |
---|
505 | e1e2v (:,:) = e1v(:,:) * e2v(:,:) ; r1_e1e2v(:,:) = 1._wp / e1e2v(:,:) |
---|
506 | |
---|
507 | e2_e1u(:,:) = e2u(:,:) / e1u(:,:) ; e1_e2v(:,:) = e1v(:,:) / e2v(:,:) |
---|
508 | !$OMP END WORKSHARE |
---|
509 | |
---|
510 | ! The following code is taken from domvvl |
---|
511 | !$OMP DO |
---|
512 | DO jk = 1, jpk |
---|
513 | DO jj = 1, jpj-1 |
---|
514 | DO ji = 1, jpi-1 |
---|
515 | !* from T- to U-point : hor. surface weighted mean |
---|
516 | e3u_n(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj) & |
---|
517 | & * ( e1e2t(ji ,jj) * ( e3t_n(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & |
---|
518 | & + e1e2t(ji+1,jj) * ( e3t_n(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) |
---|
519 | e3u_n(ji,jj,jk) = e3u_n(ji,jj,jk) + e3u_0(ji,jj,jk) |
---|
520 | |
---|
521 | !* from T- to V-point : hor. surface weighted mean |
---|
522 | e3v_n(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj) & |
---|
523 | & * ( e1e2t(ji,jj ) * ( e3t_n(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & |
---|
524 | & + e1e2t(ji,jj+1) * ( e3t_n(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) |
---|
525 | e3v_n(ji,jj,jk) = e3v_n(ji,jj,jk) + e3v_0(ji,jj,jk) |
---|
526 | END DO |
---|
527 | END DO |
---|
528 | END DO |
---|
529 | !$OMP END DO |
---|
530 | |
---|
531 | ! there needs to be an lbc_lnk call here for e3u_n and e3v_n |
---|
532 | !$OMP WORKSHARE |
---|
533 | e3u_n(jpi,:,:) = e3u_n(2,:,:) ; e3v_n(jpi,:,:) = e3v_n(2,:,:) |
---|
534 | e3u_n(:,jpj,:) = e3u_n(:,2,:) ; e3v_n(:,jpj,:) = e3v_n(:,2,:) |
---|
535 | !$OMP END WORKSHARE |
---|
536 | |
---|
537 | ! The following code is taken from domvvl: dom_vvl_init |
---|
538 | |
---|
539 | !$OMP WORKSHARE |
---|
540 | gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) |
---|
541 | gdepw_n(:,:,1) = 0.0_wp |
---|
542 | gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) |
---|
543 | !$OMP END WORKSHARE |
---|
544 | |
---|
545 | !$OMP DO PRIVATE(zcoef) |
---|
546 | DO jk = 2, jpk ! vertical sum |
---|
547 | DO jj = 1,jpj |
---|
548 | DO ji = 1,jpi |
---|
549 | ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt |
---|
550 | ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) |
---|
551 | ! ! 0.5 where jk = mikt |
---|
552 | zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) |
---|
553 | gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) |
---|
554 | gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & |
---|
555 | & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) |
---|
556 | gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) |
---|
557 | |
---|
558 | END DO |
---|
559 | END DO |
---|
560 | END DO |
---|
561 | !$OMP END DO |
---|
562 | |
---|
563 | !$OMP DO |
---|
564 | DO jj = 1, jpj |
---|
565 | DO ji = 1, jpi |
---|
566 | r1_e1u(ji,jj) = 1._wp / e1u(ji,jj) |
---|
567 | r1_e2v(ji,jj) = 1._wp / e2v(ji,jj) |
---|
568 | END DO |
---|
569 | END DO |
---|
570 | !$OMP END DO |
---|
571 | |
---|
572 | !-------------------------------------------------------------------------------------------- |
---|
573 | ! 7. Initialise fields for calculation of isopycnal diffusion |
---|
574 | |
---|
575 | !$OMP WORKSHARE |
---|
576 | uslp (:,:,:) = 0._wp ; uslpml (:,:) = 0._wp |
---|
577 | vslp (:,:,:) = 0._wp ; vslpml (:,:) = 0._wp |
---|
578 | wslpi(:,:,:) = 0._wp ; wslpiml(:,:) = 0._wp |
---|
579 | wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp |
---|
580 | |
---|
581 | akz(:,:,:) = 0._wp ; ah_wslp2(:,:,:) = 0._wp ! just to set values in the halo regions |
---|
582 | |
---|
583 | ts_b(:,:,:,:) = ts_n(:,:,:,:) ! set set ts on before time-step to the same values as the now timestep |
---|
584 | |
---|
585 | ! simplified choice of vertical diffusivity based on zdfphy |
---|
586 | avt (:,:, 1 ) = 0._wp ; avt (:,:,jpk) = 0._wp |
---|
587 | avt(:,:,2:jpkm1) = rn_avt0 |
---|
588 | !$OMP END WORKSHARE |
---|
589 | !$OMP END PARALLEL |
---|
590 | |
---|
591 | ! simplified choice of horizontal diffusivity based on ldftra: ldf_tra_init |
---|
592 | ! assumes Laplacian rather than bi-Laplacian isopycnal diffusion (inn = 1) l CASE nn_aht_ijk_t = 20 |
---|
593 | zUfac = r1_2 *rn_Ud |
---|
594 | |
---|
595 | ! from ldfc1d_c2d:ldf_c2d case TRA with knn = 1 |
---|
596 | !$OMP DO |
---|
597 | DO jj = 1, jpj |
---|
598 | DO ji = 1, jpi |
---|
599 | ahtu(ji,jj,:) = zUfac * MAX( e1u(ji,jj), e2u(ji,jj) ) |
---|
600 | ahtv(ji,jj,:) = zUfac * MAX( e1v(ji,jj), e2v(ji,jj) ) |
---|
601 | END DO |
---|
602 | END DO |
---|
603 | !$OMP END DO |
---|
604 | !$ACC END KERNELS |
---|
605 | |
---|
606 | ! taken from dommsk |
---|
607 | !$ACC KERNELS |
---|
608 | !ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 ) ! Causes a run time failure on GPU, work around below |
---|
609 | !$OMP DO |
---|
610 | DO jj = 1, jpj |
---|
611 | DO ji = 1, jpi |
---|
612 | ssmask(ji,jj)=MAXVAL(tmask(ji,jj,:), DIM=1) |
---|
613 | END DO |
---|
614 | END DO |
---|
615 | !$OMP END DO |
---|
616 | !$ACC END KERNELS |
---|
617 | |
---|
618 | !-------------------------------------------------------------------------------------------- |
---|
619 | ! 8. Calculate the pressure gradients |
---|
620 | |
---|
621 | CALL eos_init |
---|
622 | |
---|
623 | ! IF istep == 1 calculate pressure gradients without the bias pressure correction (bpc) ; if istep == 2 calculate them with the bpc |
---|
624 | |
---|
625 | WRITE(0,*)"Stepping ",maxstep,"steps" |
---|
626 | IF ( ln_zco ) THEN |
---|
627 | write(0,*) ' Calling hpg_zco' |
---|
628 | ELSE IF ( ln_zps ) THEN |
---|
629 | write(0,*) ' Calling hpg_zps' |
---|
630 | ELSE IF ( ln_sco ) THEN |
---|
631 | write(0,*) ' Calling hpg_sco' |
---|
632 | END IF |
---|
633 | IF(ln_tiled) THEN |
---|
634 | WRITE(0,*)" Tiled isopycnal diffusion",ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile |
---|
635 | ELSE |
---|
636 | WRITE(0,*)" Non-Tiled isopycnal diffusion" |
---|
637 | ENDIF |
---|
638 | WRITE(0,*) |
---|
639 | |
---|
640 | step_time = TIMER() |
---|
641 | DO istep = 1, maxstep |
---|
642 | |
---|
643 | !$ACC KERNELS |
---|
644 | IF ( istep == maxstep ) THEN |
---|
645 | !$OMP PARALLEL WORKSHARE |
---|
646 | ts_n(:,:,:,:) = ts_n(:,:,:,:) + ts_pcbias(:,:,:,:) |
---|
647 | !$OMP END PARALLEL WORKSHARE |
---|
648 | END IF |
---|
649 | |
---|
650 | !$OMP PARALLEL WORKSHARE |
---|
651 | ua(:,:,:) = 0.0_wp ; va(:,:,:) = 0.0_wp |
---|
652 | !$OMP END PARALLEL WORKSHARE |
---|
653 | |
---|
654 | !$ACC END KERNELS |
---|
655 | |
---|
656 | et = TIMER() |
---|
657 | CALL eos_insitu_pot( ts_n, tmask, rhd, rhop, gdept_n ) ! now in situ density for hpg computation; ts, tmask, gdept are IN; rhd, rhop are OUT |
---|
658 | eos_time = eos_time + (TIMER() - et) |
---|
659 | |
---|
660 | IF( ln_zps ) THEN |
---|
661 | CALL zps_hde( istep, jpts, ts_n, mbku, mbkv, e3w_n, gdept_n, tmask, umask, vmask, & ! Partial steps: before horizontal gradient |
---|
662 | & gtsu, gtsv, rhd, gru , grv, & ! of t, s, rd at the last ocean level |
---|
663 | & zti, zhi, zri, ztj, zhj, zrj ) |
---|
664 | END IF |
---|
665 | |
---|
666 | IF ( ln_zco ) THEN ! z-coordinate |
---|
667 | et = TIMER() |
---|
668 | CALL hpg_zco( grav, r1_e1u, r1_e2v, e3w_n, rhd, ua, va, zhpi ,zhpj) |
---|
669 | hpg_zco_time = hpg_zco_time + (TIMER() - et) |
---|
670 | |
---|
671 | ELSE IF ( ln_zps ) THEN ! z-coordinate plus partial steps (interpolation) |
---|
672 | et = TIMER() |
---|
673 | CALL hpg_zps( grav, r1_e1u, r1_e2v, mbku, mbkv, gru, grv, e3w_n, rhd, ua, va, zhpi ,zhpj) |
---|
674 | hpg_zps_time = hpg_zps_time + (TIMER() - et) |
---|
675 | |
---|
676 | ELSE IF ( ln_sco ) THEN |
---|
677 | et = TIMER() |
---|
678 | CALL hpg_sco( grav, r1_e1u, r1_e2v, e3w_n, rhd, gde3w_n, ln_linssh, ua, va, zhpi ,zhpj ) |
---|
679 | hpg_sco_time = hpg_sco_time + (TIMER() - et) |
---|
680 | |
---|
681 | END IF |
---|
682 | |
---|
683 | !$ACC KERNELS |
---|
684 | IF (istep == 1) THEN |
---|
685 | !$OMP PARALLEL WORKSHARE |
---|
686 | hpgi(:,:,:) = ua(:,:,:) |
---|
687 | hpgj(:,:,:) = va(:,:,:) |
---|
688 | !$OMP END PARALLEL WORKSHARE |
---|
689 | ELSE |
---|
690 | !$OMP PARALLEL WORKSHARE |
---|
691 | hbpcgi(:,:,:) = hpgi(:,:,:) - ua(:,:,:) |
---|
692 | hbpcgj(:,:,:) = hpgj(:,:,:) - va(:,:,:) |
---|
693 | !$OMP END PARALLEL WORKSHARE |
---|
694 | END IF |
---|
695 | |
---|
696 | !$OMP PARALLEL DO |
---|
697 | DO jk = 2, jpk ! vertical sum |
---|
698 | DO jj = 1,jpj |
---|
699 | DO ji = 1,jpi |
---|
700 | IF ( umask(ji,jj,jk) == 0. ) ua(ji,jj,jk) = fill |
---|
701 | IF ( vmask(ji,jj,jk) == 0. ) va(ji,jj,jk) = fill |
---|
702 | END DO |
---|
703 | END DO |
---|
704 | END DO |
---|
705 | !$OMP END PARALLEL DO |
---|
706 | !$ACC END KERNELS |
---|
707 | |
---|
708 | !$ACC KERNELS |
---|
709 | !$OMP PARALLEL WORKSHARE ! This doesn't seem to scale for OMP |
---|
710 | sum_ua=SUM(ua) |
---|
711 | sum_va=SUM(va) |
---|
712 | !$OMP END PARALLEL WORKSHARE |
---|
713 | !$ACC END KERNELS |
---|
714 | |
---|
715 | !-------------------------------------------------------------------------------------------- |
---|
716 | ! 9. Calculate isopycnal diffusion |
---|
717 | |
---|
718 | et = TIMER() |
---|
719 | |
---|
720 | CALL eos_rab_3d( ts_b, gdept_n, tmask, rab_b ) |
---|
721 | |
---|
722 | !$ACC KERNELS |
---|
723 | !$OMP PARALLEL WORKSHARE |
---|
724 | rn2b(:,:,1) = 0._wp ; rn2b(:,:,jpk) = 0._wp |
---|
725 | !$OMP END PARALLEL WORKSHARE |
---|
726 | !$ACC END KERNELS |
---|
727 | |
---|
728 | CALL bn2 ( ts_b, rab_b, gdepw_n, gdept_n, e3w_n, wmask, rn2b ) |
---|
729 | |
---|
730 | CALL zdf_mxl( rn2b, e3w_n, gdept_n, gdepw_n, avt, wmask, ssmask, mbkt, nlb10, nmln, hmld, hmlp, hmlpt ) |
---|
731 | |
---|
732 | CALL ldf_slp( ln_zps, ln_isfcav, rhd, rn2b, tmask, umask, vmask, wmask, gru, grv, & |
---|
733 | & e1t, e2t, r1_e1u, r1_e2v, e3u_n, e3v_n, e3w_n, gdept_n, gdepw_n, mbku, mbkv, & |
---|
734 | & mikt, miku, mikv, nmln, hmlp, hmlpt, risfdep, & |
---|
735 | & omlmask, uslpml, vslpml, wslpiml, wslpjml, uslp, vslp, wslpi, wslpj ) |
---|
736 | |
---|
737 | IF(.NOT.ln_tiled) THEN |
---|
738 | CALL tra_ldf_iso( ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav, jpts, kpass, r2dt, & |
---|
739 | & e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t, mbku, mbkv, miku, mikv, & |
---|
740 | & umask, vmask, wmask, e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp, & |
---|
741 | & ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts_b , ts_b, & |
---|
742 | & akz, ah_wslp2, ts_a ) |
---|
743 | ELSE |
---|
744 | !$OMP PARALLEL DO |
---|
745 | DO jsk = 1, jpkm1, jk_len_3D_tile |
---|
746 | tile % jsk = jsk |
---|
747 | tile % jsk_2 = MAX(jsk, 2) |
---|
748 | tile % jskm1 = MAX(jsk-1, 1) |
---|
749 | tile % jek = MIN(jsk + jk_len_3D_tile - 1, jpkm1) |
---|
750 | tile % jekp1 = MIN(tile % jek + 1, jpkm1) |
---|
751 | tile % jskh = tile % jsk - nn_hls |
---|
752 | tile % jekh = tile % jek + nn_hls |
---|
753 | |
---|
754 | DO jsj = nn_hls + 1, jpj - nn_hls, jj_len_3D_tile |
---|
755 | tile % jsj = jsj |
---|
756 | tile % jsjm1 = jsj-1 |
---|
757 | tile % jej = MIN(jsj + jj_len_3D_tile - 1, jpj - nn_hls) |
---|
758 | tile % jejp1 = tile % jej+1 |
---|
759 | tile % jsjh = tile % jsj - nn_hls |
---|
760 | tile % jejh = tile % jej + nn_hls |
---|
761 | DO jsi = nn_hls + 1, jpi - nn_hls, ji_len_3D_tile |
---|
762 | tile % jsi = jsi |
---|
763 | tile % jsim1 = jsi-1 |
---|
764 | tile % jei = MIN(jsi + ji_len_3D_tile - 1, jpi - nn_hls) |
---|
765 | tile % jeip1 = tile % jei+1 |
---|
766 | tile % jsih = tile % jsi - nn_hls |
---|
767 | tile % jeih = tile % jei + nn_hls |
---|
768 | |
---|
769 | CALL tra_ldf_iso_tile( tile, ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav, jpts, kpass, r2dt, & |
---|
770 | & e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t, mbku, mbkv, miku, mikv, & |
---|
771 | & umask, vmask, wmask, e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp, & |
---|
772 | & ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts_b , ts_b, ts_a, & |
---|
773 | & l_hst, l_ptr, akz, ah_wslp2, zftu, zftv ) |
---|
774 | |
---|
775 | END DO ! jis |
---|
776 | END DO ! jjs |
---|
777 | END DO ! jks |
---|
778 | !$OMP END PARALLEL DO |
---|
779 | ENDIF |
---|
780 | diffusion_time = diffusion_time + (TIMER() - et) |
---|
781 | |
---|
782 | !$ACC KERNELS |
---|
783 | !$OMP PARALLEL WORKSHARE |
---|
784 | sum_ts_a=SUM(ts_a(:,:,:,1)) |
---|
785 | !$OMP END PARALLEL WORKSHARE |
---|
786 | !$ACC END KERNELS |
---|
787 | |
---|
788 | WRITE(0,*) 'istep = ', istep |
---|
789 | WRITE(0,*) 'UA SUM ',sum_ua |
---|
790 | WRITE(0,*) 'VA SUM ',sum_va |
---|
791 | WRITE(0,*) 'TSA SUM ',sum_ts_a |
---|
792 | WRITE(0,*) 'Step time ', TIMER() - step_time |
---|
793 | WRITE(0,*) |
---|
794 | |
---|
795 | END DO ! istep |
---|
796 | |
---|
797 | step_time = TIMER() - step_time |
---|
798 | |
---|
799 | !-------------------------------------------------------------------------------------------- |
---|
800 | ! 10. Write out the pressure gradient fields |
---|
801 | |
---|
802 | write(0,*) 'writing fields: after pseudo time-steps; fields are written at tracer lats and lons for simplicity for now ' |
---|
803 | |
---|
804 | !iostat = NF90_PUT_VAR(ncid, nav_lon_varid, nav_lon) |
---|
805 | !iostat = NF90_PUT_VAR(ncid, nav_lat_varid, nav_lat) |
---|
806 | !iostat = NF90_PUT_VAR(ncid, nav_lev_varid, nav_lev) |
---|
807 | !iostat = NF90_PUT_VAR(ncid, hpgi_varid, hpgi) |
---|
808 | !iostat = NF90_PUT_VAR(ncid, hpgj_varid, hpgj) |
---|
809 | !iostat = NF90_PUT_VAR(ncid, hbpcgi_varid, hbpcgi) |
---|
810 | !iostat = NF90_PUT_VAR(ncid, hbpcgj_varid, hbpcgj) |
---|
811 | IF(ln_tiled) THEN |
---|
812 | ! !iostat = NF90_PUT_VAR(ncid, ta_varid, ts_a(:,:,:,1) ) |
---|
813 | ! !iostat = NF90_PUT_VAR(ncid, sa_varid, ts_a(:,:,:,2) ) |
---|
814 | ENDIF |
---|
815 | !IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR putting variable: ', TRIM(nf90_strerror(iostat)) |
---|
816 | ! |
---|
817 | ! Check Validation based on GPU PGI result |
---|
818 | IF ( ln_zco ) THEN ! z-coordinate |
---|
819 | WRITE(0,*)"ZCO: Difference in UA",sum_ua+409653003427.6933_wp |
---|
820 | WRITE(0,*)"ZCO: Difference in VA",sum_va+407945220514.0470_wp |
---|
821 | ELSE IF ( ln_zps ) THEN ! z-coordinate plus partial steps (interpolation) |
---|
822 | WRITE(0,*)"ZPS: Difference in UA",sum_ua+409653003427.7428_wp |
---|
823 | WRITE(0,*)"ZPS: Difference in VA",sum_va+407945220514.0659_wp |
---|
824 | ELSE IF ( ln_sco ) THEN |
---|
825 | WRITE(0,*)"SCO: Difference in UA",sum_ua+409653003427.7501_wp |
---|
826 | WRITE(0,*)"SCO: Difference in VA",sum_va+407945220514.1797_wp |
---|
827 | END IF |
---|
828 | |
---|
829 | #ifndef MANAGE |
---|
830 | !$ACC UPDATE HOST (hbpcgi, hbpcgj) |
---|
831 | #endif |
---|
832 | !-------------------------------------------------------------------------------------------- |
---|
833 | ! 9. Close the output file and stop |
---|
834 | ! |
---|
835 | iostat = NF90_PUT_VAR(ncid, hbpcgi_varid, hbpcgi) |
---|
836 | iostat = NF90_PUT_VAR(ncid, hbpcgj_varid, hbpcgj) |
---|
837 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR putting variable: ', TRIM(nf90_strerror(iostat)) |
---|
838 | iostat = NF90_CLOSE(ncid) |
---|
839 | IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing file: ', TRIM(nf90_strerror(iostat)) |
---|
840 | |
---|
841 | IF ( ln_zco ) THEN ! z-coordinate |
---|
842 | WRITE(0,*)'hpg_zco CPU time [s] is: ', hpg_zco_time |
---|
843 | ELSE IF ( ln_zps ) THEN ! z-coordinate plus partial steps (interpolation) |
---|
844 | WRITE(0,*)'hpg_zps CPU time [s] is: ', hpg_zps_time |
---|
845 | ELSE IF ( ln_sco ) THEN |
---|
846 | WRITE(0,*)'hpg_sco CPU time [s] is: ', hpg_sco_time |
---|
847 | END IF |
---|
848 | IF( ln_zps ) WRITE(0,*)'zps_hde CPU time [s] is: ', zps_hde_time |
---|
849 | WRITE(0,*)'eos_time CPU time [s] is: ', eos_time |
---|
850 | WRITE(0,*)'isopycnal diffusion CPU time [s] is: ', diffusion_time |
---|
851 | WRITE(0,*)'Total step time CPU time [s] is: ', step_time |
---|
852 | write(0,*) 'finishing' |
---|
853 | |
---|
854 | END PROGRAM main_bpc_hpg |
---|