New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
main_bpc_hpg.F90 in NEMO/branches/UKMO/BPC_miniapp/Master – NEMO

source: NEMO/branches/UKMO/BPC_miniapp/Master/main_bpc_hpg.F90 @ 10831

Last change on this file since 10831 was 10831, checked in by wayne_gaudin, 5 years ago

Ticket #2197 - added the Master directory with code and makefile

File size: 32.5 KB
Line 
1PROGRAM main_bpc_hpg
2
3!--------------------------------------------------------------------------------------------
4! 1. Declarations
5
6USE par_kind
7USE netcdf
8USE len_oce
9USE phycst
10USE in_out_manager ! I/O manager
11USE eosinsitu
12USE zpshde
13USE dynhpgs
14USE zdfmxl
15USE ldfslp
16USE traldf_iso
17USE traldf_iso_tile
18USE tiletype
19
20IMPLICIT NONE
21
22INTEGER :: ji, jj, jk, jt
23INTEGER :: jpi_in, jpj_in, jpk_in
24INTEGER :: istep,maxstep
25INTEGER :: iostat, ncid, dimid, varid, chunksize
26INTEGER :: x_dimid, y_dimid, z_dimid, i_dimid, h_dimid
27INTEGER :: x_varid, y_varid, z_varid, nav_lon_varid, nav_lat_varid, nav_lev_varid
28INTEGER :: hpgi_varid, hpgj_varid, ta_varid, sa_varid, hbpcgi_varid, hbpcgj_varid
29INTEGER :: numnam = 12  ! unit number for namelist
30INTEGER :: kpass = 1    ! for simple harmonic isopycnal diffusion
31INTEGER :: nlb10 = 7    ! w index just below 10 m depth
32
33REAL     :: fill = -99999.0   ! note that using REAL(wp) here gave errors in NF90_ENDDEF
34REAL(wp) :: r1_2 = 0.5_wp 
35REAL(wp) :: r2dt = 7200.   ! 1 hour timestep in seconds
36REAL(wp) :: zcoef, zUfac
37
38INTEGER         :: jsi, jsj, jsk   ! tile loop indices (copied into tile)
39TYPE(tile_type) :: tile 
40
41! Fields read in from files
42
43REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: in3d
44REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: in4d
45REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: nav_lon, nav_lat
46REAL(wp), DIMENSION(:),       ALLOCATABLE :: nav_lev
47REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: e1t, e1u, e1v, e2t, e2u, e2v
48REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: e3t_0, e3u_0, e3v_0, e3w_0 
49REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: tmask, umask, vmask, wmask 
50REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ts_n, ts_pcbias
51REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: e3t_n
52REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: sshn 
53REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: hpgi, hpgj, hbpcgi, hbpcgj
54
55! Fields calculated after input files have been read in
56
57REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: r1_e1u, r1_e2v, e1e2t, r1_e1e2t, e1e2u, e1e2v, r1_e1e2u, r1_e1e2v, e1_e2v, e2_e1u
58REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: e3u_n, e3v_n, e3w_n, gdept_n, gdepw_n, gde3w_n, rhd, rhop, rn2b, avt, ahtu, ahtv
59REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: gtsu, gtsv, gtui, gtvi 
60INTEGER,  DIMENSION(:,:),     ALLOCATABLE :: mbkt, mbku, mbkv, mikt, miku, mikv, nmln
61REAL(wp), DIMENSION(:,:)  ,   ALLOCATABLE :: ssmask, gru, grv, risfdep, hmld, hmlp, hmlpt, uslpml, vslpml, wslpiml, wslpjml
62REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: uslp, vslp, wslpi, wslpj, akz, ah_wslp2, omlmask
63REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ts_b, ts_a, rab_b
64REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: ua, va
65
66REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zftu, zftv    ! arrays used for diagnostics
67REAL(wp), DIMENSION(:,:)   , ALLOCATABLE  :: zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos
68REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE  :: zti, ztj             !
69REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE  :: zhpi, zhpj
70
71CHARACTER(LEN=256) :: meshfile, modelfile, pcbiasfile, outfile, out_diag 
72
73LOGICAL :: ln_linssh   ! True => linear free surface. e3t does not change with time 
74LOGICAL :: ln_zco, ln_zps, ln_sco   ! choice of hpg scheme. Only one logical should be true !
75REAL(wp):: rn_avt0, rn_Ud 
76
77LOGICAL :: ln_isfcav                ! set to false
78LOGICAL :: ln_traldf_msc, ln_traldf_blp,  ln_traldf_lap
79
80LOGICAL :: l_hst, l_ptr ! logicals for diagnostics
81
82INTEGER :: ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile
83
84LOGICAL :: ln_tiled
85
86REAL(wp) :: sum_ua,sum_va,sum_ts_a, step_time
87REAL(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
95NAMELIST/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                 
99write(0,*) 'Starting'
100
101meshfile = 'mesh_mask.nc'
102modelfile = 'restart.nc'
103pcbiasfile = 'pcbias.nc'
104outfile  = 'hpg_Jan_2019.nc'
105out_diag = 'out_diag_Jan_2019.txt'
106ln_linssh = .False.
107ln_zco  = .False. ;  ln_zps = .False. ;  ln_sco = .False. 
108
109ji_len_3D_tile = 128 ; jj_len_3D_tile = 64 ; jk_len_3D_tile = 3      ! default values for tile lengths
110
111rn_avt0 = 1.2E-5  ! values used in experiment u-bd189 (need to check what this is !!!)
112rn_ud   = 0.011   ! values used in experiment u-bd189
113
114ln_isfcav = .False. 
115
116ln_traldf_msc = .False.; ln_traldf_blp = .False.;  ln_traldf_lap = .True.
117
118l_hst = .False. ; l_ptr = .False.
119
120nlb10 = 7
121
122ln_tiled=.false.
123maxstep=10
124
125OPEN( UNIT=numnam, FILE='nam_bpc_hpg', FORM='FORMATTED', STATUS='OLD' )
126READ( numnam, nam_bpc_hpg )
127CLOSE( numnam )
128
129OPEN( UNIT=6, FILE=out_diag, FORM='FORMATTED', STATUS='UNKNOWN' )
130WRITE(6, nam_bpc_hpg ) 
131
132!--------------------------------------------------------------------------------------------
133! 3. Read in the meshmask data
134
135chunksize = 32000000
136
137write(0,*) 'Opening mesh'
138
139iostat = NF90_OPEN(TRIM(meshfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize)
140
141IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening meshfile: ', TRIM(nf90_strerror(iostat))
142
143iostat = NF90_INQ_DIMID(ncid, 'x', dimid)
144iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpi)
145iostat = NF90_INQ_DIMID(ncid, 'y', dimid)
146iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpj)
147iostat = NF90_INQ_DIMID(ncid, 'z', dimid)
148iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpk)
149IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading dimensions in meshfile: ', TRIM(nf90_strerror(iostat))
150
151write(0,*) 'jpi, jpj, jpk = ', jpi, jpj, jpk
152
153IF ( 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
157END IF
158 
159jpim1 = jpi - 1
160jpjm1 = jpj - 1
161jpkm1 = jpk - 1
162
163ALLOCATE( in3d(jpi,jpj,1), in4d(jpi,jpj,jpk,1) ) 
164
165ALLOCATE( 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
168ALLOCATE( 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
171ALLOCATE( ts_n(jpi,jpj,jpk,jpts), ts_pcbias(jpi,jpj,jpk,jpts) )
172ALLOCATE( e3t_n(jpi,jpj,jpk), sshn(jpi,jpj) ) 
173
174ALLOCATE( hpgi(jpi,jpj,jpk), hpgj(jpi,jpj,jpk), hbpcgi(jpi,jpj,jpk), hbpcgj(jpi,jpj,jpk) ) 
175
176ALLOCATE( 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
179iostat = NF90_INQ_VARID(ncid, 'nav_lon', varid)
180iostat = NF90_GET_VAR(ncid, varid, nav_lon(:,:))
181
182iostat = NF90_INQ_VARID(ncid, 'nav_lat', varid)
183iostat = NF90_GET_VAR(ncid, varid, nav_lat(:,:))
184IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading nav_lon or nav_lat: ', TRIM(nf90_strerror(iostat))
185
186iostat = NF90_INQ_VARID(ncid, 'nav_lev', varid)
187iostat = NF90_GET_VAR(ncid, varid, nav_lev(:))
188IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading nav_lev: ', TRIM(nf90_strerror(iostat))
189
190iostat = NF90_INQ_VARID(ncid, 'e1t', varid)
191iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:))
192!$OMP PARALLEL WORKSHARE
193e1t(:,:) = in3d(:,:,1)
194!$OMP END PARALLEL WORKSHARE
195
196iostat = NF90_INQ_VARID(ncid, 'e1u', varid)
197iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:))
198!$OMP PARALLEL WORKSHARE
199e1u(:,:) = in3d(:,:,1)
200!$OMP END PARALLEL WORKSHARE
201
202iostat = NF90_INQ_VARID(ncid, 'e1v', varid)
203iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:))
204!$OMP PARALLEL WORKSHARE
205e1v(:,:) = in3d(:,:,1)
206!$OMP END PARALLEL WORKSHARE
207
208iostat = NF90_INQ_VARID(ncid, 'e2t', varid)
209iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:))
210!$OMP PARALLEL WORKSHARE
211e2t(:,:) = in3d(:,:,1)
212!$OMP END PARALLEL WORKSHARE
213
214iostat = NF90_INQ_VARID(ncid, 'e2u', varid)
215iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:))
216!$OMP PARALLEL WORKSHARE
217e2u(:,:) = in3d(:,:,1)
218!$OMP END PARALLEL WORKSHARE
219
220iostat = NF90_INQ_VARID(ncid, 'e2v', varid)
221iostat = NF90_GET_VAR(ncid, varid, in3d(:,:,:))
222!$OMP PARALLEL WORKSHARE
223e2v(:,:) = in3d(:,:,1)
224!$OMP END PARALLEL WORKSHARE
225
226IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e1, e2 fields: ', TRIM(nf90_strerror(iostat))
227
228iostat = NF90_INQ_VARID(ncid, 'e3t_0', varid)
229iostat = NF90_GET_VAR(ncid, varid, e3t_0(:,:,:))
230
231iostat = NF90_INQ_VARID(ncid, 'e3u_0', varid)
232iostat = NF90_GET_VAR(ncid, varid, e3u_0(:,:,:))
233
234iostat = NF90_INQ_VARID(ncid, 'e3v_0', varid)
235iostat = NF90_GET_VAR(ncid, varid, e3v_0(:,:,:))
236
237iostat = NF90_INQ_VARID(ncid, 'e3w_0', varid)
238iostat = NF90_GET_VAR(ncid, varid, e3w_0(:,:,:))
239
240IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading e3 fields: ', TRIM(nf90_strerror(iostat))
241
242iostat = NF90_INQ_VARID(ncid, 'tmask', varid)
243iostat = NF90_GET_VAR(ncid, varid, tmask(:,:,:))
244
245iostat = NF90_INQ_VARID(ncid, 'umask', varid)
246iostat = NF90_GET_VAR(ncid, varid, umask(:,:,:))
247
248iostat = NF90_INQ_VARID(ncid, 'vmask', varid)
249iostat = NF90_GET_VAR(ncid, varid, vmask(:,:,:))
250
251IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading mask fields: ', TRIM(nf90_strerror(iostat))
252
253iostat = NF90_CLOSE(ncid)
254IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing mesh file: ', TRIM(nf90_strerror(iostat))
255
256!--------------------------------------------------------------------------------------------
257! 4. Read in the models fields required
258
259write(0,*) 'Opening model fields'
260
261iostat = NF90_OPEN(TRIM(modelfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize)
262IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening main fields file: ', TRIM(nf90_strerror(iostat))
263
264iostat = NF90_INQ_DIMID(ncid, 'x', dimid)
265iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpi)
266iostat = NF90_INQ_DIMID(ncid, 'y', dimid)
267iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpj)
268iostat = NF90_INQ_DIMID(ncid, 'z', dimid)
269iostat = NF90_INQUIRE_DIMENSION(ncid, dimid, LEN=jpk)
270IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading main field dimensions: ', TRIM(nf90_strerror(iostat))
271
272write(0,*) 'jpi, jpj, jpk = ', jpi, jpj, jpk
273
274IF ( 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
278END IF
279
280iostat = NF90_INQ_VARID(ncid, 'tn', varid)
281iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:))
282!$OMP PARALLEL WORKSHARE
283ts_n(:,:,:,1) = in4d(:,:,:,1)
284!$OMP END PARALLEL WORKSHARE
285IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_n(:,:,:,1) field: ', TRIM(nf90_strerror(iostat))
286
287iostat = NF90_INQ_VARID(ncid, 'sn', varid)
288iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:))
289!$OMP PARALLEL WORKSHARE
290ts_n(:,:,:,2) = in4d(:,:,:,1)
291!$OMP END PARALLEL WORKSHARE
292IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_n(:,:,:,2) field: ', TRIM(nf90_strerror(iostat))
293
294IF ( .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))
301ELSE
302!$OMP PARALLEL WORKSHARE
303   e3t_n(:,:,:) = e3t_0(:,:,:)
304!$OMP END PARALLEL WORKSHARE
305END IF
306
307IF ( 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))
314ELSE 
315!$OMP PARALLEL WORKSHARE
316   sshn(:,:) = 0._wp 
317!$OMP END PARALLEL WORKSHARE
318END IF
319
320iostat = NF90_CLOSE(ncid)
321IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing main fields (file: ', TRIM(nf90_strerror(iostat))
322
323write(0,*) 'Opening pcbias fields'
324
325iostat = NF90_OPEN(TRIM(pcbiasfile), MODE=nf90_nowrite, NCID=ncid, CHUNKSIZE=chunksize)
326IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR opening pcbias (file: ', TRIM(nf90_strerror(iostat))
327
328iostat = NF90_INQ_VARID(ncid, 'tbias_p', varid)
329iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:))
330!$OMP PARALLEL WORKSHARE
331ts_pcbias(:,:,:,1) = in4d(:,:,:,1)
332!$OMP END PARALLEL WORKSHARE
333IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_pcbias(:,:,:,1) field: ', TRIM(nf90_strerror(iostat))
334
335iostat = NF90_INQ_VARID(ncid, 'sbias_p', varid)
336iostat = NF90_GET_VAR(ncid, varid, in4d(:,:,:,:))
337!$OMP PARALLEL WORKSHARE
338ts_pcbias(:,:,:,2) = in4d(:,:,:,1)
339!$OMP END PARALLEL WORKSHARE
340IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR reading ts_pcbias(:,:,:,2) field: ', TRIM(nf90_strerror(iostat))
341
342iostat = NF90_CLOSE(ncid)
343IF (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
347iostat = NF90_CREATE(TRIM(outfile), NF90_CLOBBER, ncid)
348IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR creating file: ', TRIM(nf90_strerror(iostat))
349
350iostat = NF90_DEF_DIM(ncid, 'x', jpi, x_dimid)
351IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat))
352iostat = NF90_DEF_DIM(ncid, 'y', jpj, y_dimid)
353IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat))
354iostat = NF90_DEF_DIM(ncid, 'z', jpk, z_dimid)
355IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining dimension: ', TRIM(nf90_strerror(iostat))
356
357iostat = NF90_DEF_VAR(ncid, 'nav_lon', NF90_FLOAT, (/ x_dimid,y_dimid /), nav_lon_varid)
358IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat))
359iostat = NF90_DEF_VAR(ncid, 'nav_lat', NF90_FLOAT, (/ x_dimid,y_dimid /), nav_lat_varid)
360IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat))
361iostat = NF90_DEF_VAR(ncid, 'nav_lev', NF90_FLOAT, (/ z_dimid /), nav_lev_varid)
362IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining nav_lon: ', TRIM(nf90_strerror(iostat))
363
364iostat = NF90_DEF_VAR(ncid, 'hpgi', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hpgi_varid)
365IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hpgi: ', TRIM(nf90_strerror(iostat))
366iostat = NF90_DEF_VAR(ncid, 'hpgj', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hpgj_varid)
367IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hpgj: ', TRIM(nf90_strerror(iostat))
368iostat = NF90_DEF_VAR(ncid, 'hbpcgi', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hbpcgi_varid)
369IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hbpcgi: ', TRIM(nf90_strerror(iostat))
370iostat = NF90_DEF_VAR(ncid, 'hbpcgj', NF90_FLOAT, (/ x_dimid,y_dimid,z_dimid /), hbpcgj_varid)
371IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining hbpcgj: ', TRIM(nf90_strerror(iostat))
372
373iostat = NF90_PUT_ATT(ncid, hpgi_varid, '_FillValue', fill)
374IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining fill hpgi_varid : ', TRIM(nf90_strerror(iostat))
375iostat = NF90_PUT_ATT(ncid, hpgj_varid, '_FillValue', fill)
376iostat = NF90_PUT_ATT(ncid, hbpcgi_varid, '_FillValue', fill)
377iostat = NF90_PUT_ATT(ncid, hbpcgj_varid, '_FillValue', fill)
378IF(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)
383ENDIF
384IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR defining fill hbpcgj_varid : ', TRIM(nf90_strerror(iostat))
385
386iostat = NF90_ENDDEF(ncid)
387IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR leaving define mode: ', TRIM(nf90_strerror(iostat))
388!set timers
389hpg_zco_time = 0.
390hpg_zps_time = 0.
391hpg_sco_time = 0.
392zps_hde_time = 0.
393eos_time     = 0.
394eos2d_time   = 0.
395diffusion_time = 0.
396step_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
409ALLOCATE( 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
412ALLOCATE( e3u_n(jpi,jpj,jpk), e3v_n(jpi,jpj,jpk), e3w_n(jpi,jpj,jpk) )
413
414ALLOCATE( gdept_n(jpi,jpj,jpk), gdepw_n(jpi,jpj,jpk), gde3w_n(jpi,jpj,jpk) )
415
416ALLOCATE( 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
418ALLOCATE( gtsu(jpi,jpj, jpts), gtsv(jpi,jpj, jpts), gtui(jpi,jpj,jpts),  gtvi(jpi,jpj,jpts) ) 
419
420ALLOCATE( mbkt(jpi,jpj), mbku(jpi,jpj), mbkv(jpi,jpj), mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), nmln(jpi,jpj) )
421
422ALLOCATE( 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
426ALLOCATE( uslp(jpi,jpj,jpk), vslp(jpi,jpj,jpk), wslpi(jpi,jpj,jpk), wslpj(jpi,jpj,jpk) )   
427
428ALLOCATE(  akz(jpi,jpj,jpk), ah_wslp2(jpi,jpj,jpk), omlmask(jpi,jpj,jpk) ) 
429
430ALLOCATE( ts_b(jpi,jpj,jpk,jpts), ts_a(jpi,jpj,jpk,jpts), rab_b(jpi,jpj,jpk,jpts) )
431
432ALLOCATE( ua(jpi,jpj,jpk), va(jpi,jpj,jpk) )
433
434IF(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
445DO 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
453END DO 
454!$OMP END DO
455
456!$OMP DO
457DO 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
463END DO
464!$OMP END DO
465
466! set ice-shelf fields for no ice !
467mikt(:,:) = 1 ; miku(:,:) = 1 ; mikv(:,:) = 1  ; risfdep(:,:) = 0._wp 
468gtui(:,:,:) = 0._wp ; gtvi(:,:,:) = 0._wp
469 
470! The following code is taken from dommsk: dom_msk
471!$OMP WORKSHARE
472wmask (:,:,1) = tmask(:,:,1)
473!$OMP END WORKSHARE
474!$OMP DO
475DO jk = 2, jpk                   ! interior values
476   wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
477END DO 
478!$OMP END DO
479
480! The following code is taken from domvvl: dom_vvl_interpol for case 'W' interpolation
481
482IF ( ln_linssh ) THEN 
483!$OMP WORKSHARE
484  e3w_n(:,:,:) = e3w_0(:,:,:)
485!$OMP END WORKSHARE
486
487ELSE 
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
499END IF 
500
501! The following code is taken from domhgr: SUBROUTINE dom_hgr
502!$OMP WORKSHARE
503e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
504e1e2u (:,:) = e1u(:,:) * e2u(:,:)   ;   r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)
505e1e2v (:,:) = e1v(:,:) * e2v(:,:)   ;   r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
506
507e2_e1u(:,:) = e2u(:,:) / e1u(:,:)   ;   e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
508!$OMP END WORKSHARE
509
510! The following code is taken from domvvl
511!$OMP DO
512DO 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
528END 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
533e3u_n(jpi,:,:) = e3u_n(2,:,:) ; e3v_n(jpi,:,:) = e3v_n(2,:,:)
534e3u_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
540gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
541gdepw_n(:,:,1) = 0.0_wp
542gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
543!$OMP END WORKSHARE
544
545!$OMP DO PRIVATE(zcoef)
546DO 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
560END DO
561!$OMP END DO
562
563!$OMP DO
564DO 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
569END DO 
570!$OMP END DO
571
572!--------------------------------------------------------------------------------------------
573! 7. Initialise fields for calculation of isopycnal diffusion   
574
575!$OMP WORKSHARE
576uslp (:,:,:) = 0._wp   ;   uslpml (:,:) = 0._wp     
577vslp (:,:,:) = 0._wp   ;   vslpml (:,:) = 0._wp
578wslpi(:,:,:) = 0._wp   ;   wslpiml(:,:) = 0._wp
579wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp
580
581akz(:,:,:) = 0._wp     ;   ah_wslp2(:,:,:) = 0._wp    ! just to set values in the halo regions
582
583ts_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
586avt  (:,:, 1 ) = 0._wp  ; avt  (:,:,jpk) = 0._wp
587avt(:,:,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
593zUfac = r1_2 *rn_Ud
594
595! from ldfc1d_c2d:ldf_c2d case  TRA with knn = 1
596!$OMP DO
597DO 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
602END 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
610DO jj = 1, jpj
611  DO ji = 1, jpi
612    ssmask(ji,jj)=MAXVAL(tmask(ji,jj,:), DIM=1)
613  END DO
614END DO 
615!$OMP END DO
616!$ACC END KERNELS
617
618!--------------------------------------------------------------------------------------------
619! 8. Calculate the pressure gradients 
620
621CALL 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
625WRITE(0,*)"Stepping ",maxstep,"steps"
626IF ( ln_zco ) THEN
627   write(0,*) '  Calling hpg_zco'
628ELSE IF ( ln_zps ) THEN
629   write(0,*) '  Calling hpg_zps'
630ELSE IF ( ln_sco ) THEN
631   write(0,*) '  Calling hpg_sco'
632END IF
633IF(ln_tiled) THEN
634  WRITE(0,*)"  Tiled isopycnal diffusion",ji_len_3D_tile, jj_len_3D_tile, jk_len_3D_tile
635ELSE
636  WRITE(0,*)"  Non-Tiled isopycnal diffusion"
637ENDIF
638WRITE(0,*)
639
640step_time = TIMER()
641DO 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
795END DO ! istep
796
797step_time = TIMER() - step_time
798
799!--------------------------------------------------------------------------------------------
800! 10. Write out the pressure gradient fields
801     
802write(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)
811IF(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) )
814ENDIF
815!IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR putting variable: ', TRIM(nf90_strerror(iostat))
816!
817! Check Validation based on GPU PGI result
818IF ( 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
821ELSE 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
824ELSE 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
827END IF 
828
829#ifndef MANAGE
830!$ACC UPDATE HOST (hbpcgi, hbpcgj)
831#endif
832!--------------------------------------------------------------------------------------------
833! 9. Close the output file and stop
834!
835iostat = NF90_PUT_VAR(ncid, hbpcgi_varid, hbpcgi)
836iostat = NF90_PUT_VAR(ncid, hbpcgj_varid, hbpcgj)
837IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR putting variable: ', TRIM(nf90_strerror(iostat))
838iostat = NF90_CLOSE(ncid)
839IF (iostat /= nf90_noerr) WRITE(numout,*) 'ERROR closing file: ', TRIM(nf90_strerror(iostat))
840
841IF ( ln_zco ) THEN    ! z-coordinate
842  WRITE(0,*)'hpg_zco  CPU time [s] is:             ', hpg_zco_time
843ELSE IF ( ln_zps ) THEN  ! z-coordinate plus partial steps (interpolation)
844  WRITE(0,*)'hpg_zps  CPU time [s] is:             ', hpg_zps_time
845ELSE IF ( ln_sco ) THEN
846  WRITE(0,*)'hpg_sco  CPU time [s] is:             ', hpg_sco_time
847END IF
848IF( ln_zps ) WRITE(0,*)'zps_hde  CPU time [s] is:             ', zps_hde_time
849WRITE(0,*)'eos_time CPU time [s] is:             ', eos_time
850WRITE(0,*)'isopycnal diffusion  CPU time [s] is: ', diffusion_time
851WRITE(0,*)'Total step time      CPU time [s] is: ', step_time
852write(0,*) 'finishing'
853
854END PROGRAM main_bpc_hpg
Note: See TracBrowser for help on using the repository browser.