Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diahsb_gpu.h90 @ 15063

Last change on this file since 15063 was 14815, checked in by mcastril, 3 years ago
[2021/HPC-11_mcastril_HPDAonline_DiagGPU] Final update for diahsb and diahsb_gpu
File size: 32.7 KB

Line
1	MODULE diahsb
2	!!======================================================================
3	!! * MODULE diahsb *
4	!! Ocean diagnostics: Heat, salt and volume budgets
5	!!======================================================================
6	!! History : 3.3 ! 2010-09 (M. Leclair) Original code
7	!! ! 2012-10 (C. Rousset) add iom_put
8	!!----------------------------------------------------------------------
9
10	!!----------------------------------------------------------------------
11	!! dia_hsb : Diagnose the conservation of ocean heat and salt contents, and volume
12	!! dia_hsb_rst : Read or write DIA file in restart file
13	!! dia_hsb_init : Initialization of the conservation diagnostic
14	!!----------------------------------------------------------------------
15	USE oce ! ocean dynamics and tracers
16	USE dom_oce ! ocean space and time domain
17	USE phycst ! physical constants
18	USE sbc_oce ! surface thermohaline fluxes
19	USE isf_oce ! ice shelf fluxes
20	USE sbcrnf ! river runoff
21	USE domvvl ! vertical scale factors
22	USE traqsr ! penetrative solar radiation
23	USE trabbc ! bottom boundary condition
24	USE trabbc ! bottom boundary condition
25	USE restart ! ocean restart
26	USE bdy_oce , ONLY : ln_bdy
27	!
28	USE iom ! I/O manager
29	USE in_out_manager ! I/O manager
30	USE gpu_manager ! GPU manager
31	USE cudafor ! CUDA toolkit libs
32	USE cuda_fortran ! CUDA routines
33	!USE nvtx ! CUDA profiling/DEGUG tools
34	USE lib_fortran ! glob_sum
35	USE lib_mpp ! distributed memory computing library
36	USE timing ! preformance summary
37
38	IMPLICIT NONE
39	PRIVATE
40
41	PUBLIC dia_hsb ! routine called by step.F90
42	PUBLIC dia_hsb_init ! routine called by nemogcm.F90
43
44	LOGICAL, PUBLIC :: ln_diahsb !: check the heat and salt budgets
45
46	REAL(wp) :: surf_tot ! ocean surface
47	REAL(wp) , DIMENSION(2), SAVE :: frc_t, frc_s, frc_v ! global forcing trends
48	REAL(wp) , DIMENSION(2), SAVE :: frc_wn_t, frc_wn_s ! global forcing trends
49	!
50	REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf
51	REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf_ini , ssh_ini !
52	REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: ssh_hc_loc_ini, ssh_sc_loc_ini !
53	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, PINNED :: hc_loc_ini, sc_loc_ini !
54	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3t_ini !
55	REAL(wp), DIMENSION(:) , ALLOCATABLE, PINNED, SAVE :: h_ztmpv, h_ztmph, h_ztmps, h_ztmp !
56	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmask_ini
57
58
59	!Device data associate to PUBLIC arrays
60	REAL(8), DIMENSION(:,:,:,:) , ALLOCATABLE, DEVICE :: d_e3t !
61	REAL(8), DIMENSION(:,:,:) , ALLOCATABLE, DEVICE :: d_tmask !
62	REAL(8), DIMENSION(:,:) , ALLOCATABLE, DEVICE :: d_tmask_h !
63	REAL(8), DIMENSION(:,:,:) , ALLOCATABLE, DEVICE :: d_tmask_ini !
64	REAL(8), DIMENSION(:,:,:,:,:), ALLOCATABLE, DEVICE :: d_ts !
65	!Device data associate to LOCAL/DEVICE arrays
66	REAL(8), DEVICE , DIMENSION(:,:) , ALLOCATABLE :: d_surf !
67	REAL(8), DEVICE , DIMENSION(:,:) , ALLOCATABLE :: d_surf_ini !
68	REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_hc_loc_ini !
69	REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_sc_loc_ini !
70	REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_e3t_ini !
71	REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_zwrkv, d_zwrkh, d_zwrks, d_zwrk ! 3D GPU workspace
72	REAL(8), DEVICE :: ztmpv, ztmph, ztmps, ztmp ! Device Reduction
73	!
74	INTEGER :: globsize ! 3D workspace size
75	type(dim3) :: dimGrid, dimBlock ! cuda parameters
76	INTEGER, parameter :: nstreams = 3 ! Streams Number
77	INTEGER(kind=cuda_stream_kind) :: stream(nstreams), str ! Stream ID
78	!DEBUG
79	!REAL(8) , save , DIMENSION(:,:,:) , ALLOCATABLE :: prev_3d
80	!REAL(8) :: accum
81
82
83
84
85	!! * Substitutions
86	# include "domzgr_substitute.h90"
87	!!----------------------------------------------------------------------
88	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
89	!! $Id$
90	!! Software governed by the CeCILL license (see ./LICENSE)
91	!!----------------------------------------------------------------------
92	CONTAINS
93
94	SUBROUTINE dia_hsb( kt, Kbb, Kmm )
95	!!---------------------------------------------------------------------------
96	!! * ROUTINE dia_hsb *
97	!!
98	!! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
99	!!
100	!! ** Method : - Compute the deviation of heat content, salt content and volume
101	!! at the current time step from their values at nit000
102	!! - Compute the contribution of forcing and remove it from these deviations
103	!!
104	!!---------------------------------------------------------------------------
105	INTEGER, INTENT(in) :: kt ! ocean time-step index
106	INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices
107	!
108	INTEGER, VALUE :: ji, jj, jk, kts ! dummy loop indice
109	INTEGER, VALUE :: localsize ! jpi * jpj * jpk
110	INTEGER :: istat ! CUDA error check
111	COMPLEX :: ctmp ! dummy complex number
112	INTEGER(kind=cuda_stream_kind) :: str ! dummy kernel stream
113	INTEGER :: tile_n, tile_b ! tile indexe. _n now, _b before
114	REAL(wp) , DIMENSION(2), SAVE :: zdiff_hc1, zdiff_sc1 ! heat and salt content variations
115	REAL(wp) , DIMENSION(2), SAVE :: zdiff_hc, zdiff_sc ! - - - -
116	REAL(wp) , DIMENSION(2), SAVE :: zdiff_v2 ! volume variation
117	REAL(wp) , DIMENSION(2), SAVE :: zdiff_v1 ! volume variation
118	REAL(wp) , DIMENSION(2), SAVE :: zerr_hc1, zerr_sc1 ! heat and salt content misfit
119	REAL(wp) , DIMENSION(2), SAVE :: zvol_tot ! volume
120	REAL(wp) , DIMENSION(2), SAVE :: z_frc_trd_t, z_frc_trd_s ! - -
121	REAL(wp) , DIMENSION(2), SAVE :: z_frc_trd_v ! - -
122	REAL(wp) , DIMENSION(2), SAVE :: z_wn_trd_t, z_wn_trd_s ! - -
123	REAL(wp) , DIMENSION(2), SAVE :: z_ssh_hc, z_ssh_sc ! - -
124	# 147 "diahsb_new.F90"
125	REAL(wp), DIMENSION(jpi,jpj) :: z2d0, z2d1 ! 2D workspace
126	REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwrk ! 3D workspace
127	!!---------------------------------------------------------------------------
128	IF( ln_timing ) CALL timing_start('dia_hsb')
129
130	localsize = jpi * jpj * jpk
131	kts = kt
132	IF (kts == 1) THEN
133	tile_n = 1
134	tile_b = 1
135	ELSE
136	IF( MOD(kts,2) == 0) THEN
137	tile_n = 2
138	tile_b = 1
139	ELSE IF( MOD(kts,2) == 1 ) THEN
140	tile_n = 1
141	tile_b = 2
142	END IF
143	END IF
144
145	!
146	ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
147	ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
148	! ------------------------- !
149	! 1 - Trends due to forcing !
150	! ------------------------- !
151	z_frc_trd_v(tile_n) = r1_rho0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )! volume fluxes
152	z_frc_trd_t(tile_n) = glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) ) ! heat fluxes
153	z_frc_trd_s(tile_n) = glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) ) ! salt fluxes
154	! ! Add runoff heat & salt input
155	IF( ln_rnf ) z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) + glob_sum( 'diahsb', rnf_tsc(:,:,jp_tem) * surf(:,:) )
156	IF( ln_rnf_sal) z_frc_trd_s(tile_n) = z_frc_trd_s(tile_n) + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) ) ! Add ice shelf heat & salt input
157	! ! Add ice shelf heat & salt input
158	IF( ln_isf ) z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) &
159	& + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ) ! Add penetrative solar radiation
160	! ! Add penetrative solar radiation
161	IF( ln_traqsr ) z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) + r1_rho0_rcp * glob_sum( 'diahsb', qsr (:,:) * surf(:,:) ) ! Add geothermal heat flux
162	! ! Add geothermal heat flux
163	IF( ln_trabbc ) z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) + glob_sum( 'diahsb', qgh_trd0(:,:) * surf(:,:) )
164	!
165	IF( ln_linssh ) THEN
166	IF( ln_isfcav ) THEN
167	DO ji=1,jpi
168	DO jj=1,jpj
169	z2d0(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
170	z2d1(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
171	END DO
172	END DO
173	ELSE
174	z2d0(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
175	z2d1(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
176	END IF
177	z_wn_trd_t(tile_n) = - glob_sum( 'diahsb', z2d0 )
178	z_wn_trd_s(tile_n) = - glob_sum( 'diahsb', z2d1 )
179	ENDIF
180
181	IF (kts>1) THEN
182	frc_v(tile_n) = frc_v(tile_b)
183	frc_t(tile_n) = frc_t(tile_b)
184	frc_s(tile_n) = frc_s(tile_b)
185	frc_wn_t(tile_n) = frc_wn_t(tile_b)
186	frc_wn_s(tile_n) = frc_wn_s(tile_b)
187	END IF
188	frc_v(tile_n) = frc_v(tile_n) + z_frc_trd_v(tile_n) * rn_Dt
189	frc_t(tile_n) = frc_t(tile_n) + z_frc_trd_t(tile_n) * rn_Dt
190	frc_s(tile_n) = frc_s(tile_n) + z_frc_trd_s(tile_n) * rn_Dt
191	! ! Advection flux through fixed surface (z=0)
192	IF( ln_linssh ) THEN
193	frc_wn_t(tile_n) = frc_wn_t(tile_n) + z_wn_trd_t(tile_n) * rn_Dt
194	frc_wn_s(tile_n) = frc_wn_s(tile_n) + z_wn_trd_s(tile_n) * rn_Dt
195	ENDIF
196
197	! ------------------------ !
198	! 2 - Content variations !
199	! ------------------------ !
200	! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
201
202	! ! volume variation (calculated with ssh)
203
204	zdiff_v1(tile_n) = glob_sum_full( 'diahsb', surf(:,:)ssh(:,:,Kmm) - surf_ini(:,:)ssh_ini(:,:) )
205
206	! ! heat & salt content variation (associated with ssh)
207	IF( ln_linssh ) THEN ! linear free surface case
208	IF( ln_isfcav ) THEN ! ISF case
209	DO ji = 1, jpi
210	DO jj = 1, jpj
211	z2d0(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
212	z2d1(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
213	END DO
214	END DO
215	ELSE ! no under ice-shelf seas
216	z2d0(:,:) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )
217	z2d1(:,:) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )
218	END IF
219	z_ssh_hc(tile_n) = glob_sum_full( 'diahsb', z2d0 )
220	z_ssh_sc(tile_n) = glob_sum_full( 'diahsb', z2d1 )
221	ENDIF
222
223	str = stream(tile_n)
224	istat = 0
225	istat = cudaMemcpyAsync( d_e3t, e3t , jpijpjjpk*jpt , str ) + istat
226	istat = cudaMemcpyAsync(d_hc_loc_ini, hc_loc_ini, jpijpjjpk , str ) + istat
227	istat = cudaMemcpyAsync(d_sc_loc_ini, sc_loc_ini, jpijpjjpk , str ) + istat
228	istat = cudaMemcpyAsync( d_ts, ts , jpijpjjpk2jpt, str ) + istat
229	IF( istat /= 0 ) THEN
230	CALL ctl_stop( 'dia_hsb: unable to async GPU copy H2D' ) ; RETURN
231	ENDIF
232	dimBlock = dim3(4,4,4)
233	dimGrid = dim3( ceiling( real( jpi ) / dimBlock%x ) , ceiling( real( jpj ) / dimBlock%y ) , &
234	ceiling( real( jpkm1 ) / dimBlock%z ) )
235	!
236	CALL dia_hsb_kernel<<<dimGrid, dimBlock, 0, str>>> (d_surf , d_e3t, d_surf_ini, d_e3t_ini, &
237	& d_ts, d_hc_loc_ini, d_sc_loc_ini, d_tmask, d_tmask_ini, d_zwrkv, d_zwrkh, d_zwrks, d_zwrk, jpi, jpj, jpk, jpt, Kmm)
238
239
240	CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrkv , d_tmask_h , jpi, jpj, jpk)
241	CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrkh , d_tmask_h , jpi, jpj, jpk)
242	CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrks , d_tmask_h , jpi, jpj, jpk)
243	CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrk , d_tmask_h , jpi, jpj, jpk)
244
245	ztmpv = 0.e0
246	ztmph = 0.e0
247	ztmps = 0.e0
248	ztmp = 0.e0
249
250	istat = 0
251	!$cuf kernel do <<< , , stream=str >>>
252	do ji = 1, localsize
253	ztmpv = ztmpv + d_zwrkv(ji)
254	end do
255	istat = cudaMemcpyAsync( h_ztmpv(tile_n) , ztmpv , 1 , str ) + istat
256	!$cuf kernel do <<< , , stream=str >>>
257	do ji = 1, localsize
258	ztmph = ztmph + d_zwrkh(ji)
259	end do
260	istat = cudaMemcpyAsync( h_ztmph(tile_n) , ztmph , 1 , str ) + istat
261	!$cuf kernel do <<< , , stream=str >>>
262	do ji = 1, localsize
263	ztmps = ztmps + d_zwrks(ji)
264	end do
265	istat = cudaMemcpyAsync( h_ztmps(tile_n) , ztmps , 1 , str ) + istat
266	!$cuf kernel do <<< , , stream=str >>>
267	do ji = 1, localsize
268	ztmp = ztmp + d_zwrk(ji)
269	end do
270	istat = cudaMemcpyAsync( h_ztmp (tile_n) , ztmp , 1 , str ) + istat
271	!
272	IF( istat /= 0 ) THEN
273	CALL ctl_stop( 'dia_hsb: unable to async GPU copy D2H' ) ; RETURN
274	ENDIF
275	!
276	istat = cudaStreamSynchronize(stream(tile_b))
277	!
278	IF( istat /= 0 ) THEN
279	CALL ctl_stop( 'dia_hsb: unable to stream synchronize' ) ; RETURN
280	ENDIF
281	!
282	ctmp = CMPLX( h_ztmpv(tile_b) , 0.e0, 8 )
283	CALL mpp_sum('diahsb', ctmp )
284	zdiff_v2(tile_b) = REAL( ctmp, 8 )
285
286	ctmp = CMPLX( h_ztmph(tile_b) , 0.e0, 8 )
287	CALL mpp_sum('diahsb', ctmp )
288	zdiff_hc(tile_b) = REAL( ctmp, 8 )
289
290	ctmp = CMPLX( h_ztmps(tile_b) , 0.e0, 8 )
291	CALL mpp_sum('diahsb', ctmp )
292	zdiff_sc(tile_b) = REAL( ctmp, 8 )
293
294	ctmp = CMPLX( h_ztmp(tile_b) , 0.e0, 8 )
295	CALL mpp_sum('diahsb', ctmp )
296	zvol_tot(tile_b) = REAL( ctmp, 8 )
297
298	IF ( kt == nitend ) THEN
299	!
300	istat = cudaStreamSynchronize(stream(tile_n))
301	IF( istat /= 0 ) THEN
302	CALL ctl_stop( 'dia_hsb: unable to stream synchronize' ) ; RETURN
303	ENDIF
304	!
305	ctmp = CMPLX( h_ztmpv(tile_n) , 0.e0, 8 )
306	CALL mpp_sum('diahsb', ctmp )
307	zdiff_v2(tile_n) = REAL( ctmp, 8 )
308
309	ctmp = CMPLX( h_ztmph(tile_n) , 0.e0, 8 )
310	CALL mpp_sum('diahsb', ctmp )
311	zdiff_hc(tile_n) = REAL( ctmp, 8 )
312
313	ctmp = CMPLX( h_ztmps(tile_n) , 0.e0, 8 )
314	CALL mpp_sum('diahsb', ctmp )
315	zdiff_sc(tile_n) = REAL( ctmp, 8 )
316
317	ctmp = CMPLX( h_ztmp(tile_n) , 0.e0, 8 )
318	CALL mpp_sum('diahsb', ctmp )
319	zvol_tot(tile_n) = REAL( ctmp, 8 )
320	ENDIF
321	! ------------------------ !
322	! 3 - Drifts !
323	! ------------------------ !
324
325	IF ( kt > 1 ) THEN
326	kts = kts - 1
327	zdiff_v1(tile_b) = zdiff_v1(tile_b) - frc_v(tile_b)
328	IF( .NOT.ln_linssh ) zdiff_v2(tile_b) = zdiff_v2(tile_b) - frc_v(tile_b)
329	zdiff_hc(tile_b) = zdiff_hc(tile_b) - frc_t(tile_b)
330	zdiff_sc(tile_b) = zdiff_sc(tile_b) - frc_s(tile_b)
331	IF( ln_linssh ) THEN
332	zdiff_hc1(tile_b) = zdiff_hc (tile_b) + z_ssh_hc(tile_b)
333	zdiff_sc1(tile_b) = zdiff_sc (tile_b) + z_ssh_sc(tile_b)
334	zerr_hc1 (tile_b) = z_ssh_hc(tile_b) - frc_wn_t(tile_b)
335	zerr_sc1 (tile_b) = z_ssh_sc(tile_b) - frc_wn_s(tile_b)
336	ENDIF
337	!!gm to be added ?
338	! IF( ln_linssh ) THEN ! fixed volume, add the ssh contribution
339	! zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * sshn(:,:) )
340	! ENDIF
341	!!gm end
342
343	CALL iom_put( 'bgfrcvol' , frc_v(tile_b) * 1.e-9 ) ! vol - surface forcing (km3)
344	CALL iom_put( 'bgfrctem' , frc_t(tile_b) * rho0 * rcp * 1.e-20 ) ! hc - surface forcing (1.e20 J)
345	CALL iom_put( 'bgfrchfx' , frc_t(tile_b) * rho0 * rcp / & ! hc - surface forcing (W/m2)
346	& ( surf_tot * kts * rn_Dt ) )
347	CALL iom_put( 'bgfrcsal' , frc_s(tile_b) * 1.e-9 ) ! sc - surface forcing (psu*km3)
348	IF( .NOT. ln_linssh ) THEN
349	CALL iom_put( 'bgtemper' , zdiff_hc(tile_b) / zvol_tot(tile_b) ) ! Temperature drift (C)
350	CALL iom_put( 'bgsaline' , zdiff_sc(tile_b) / zvol_tot(tile_b) ) ! Salinity drift (PSU)
351	CALL iom_put( 'bgheatco' , zdiff_hc(tile_b) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
352	CALL iom_put( 'bgheatfx' , zdiff_hc(tile_b) * rho0 * rcp / & ! Heat flux drift (W/m2)
353	& ( surf_tot * kts * rn_Dt ) )
354	CALL iom_put( 'bgsaltco' , zdiff_sc(tile_b) * 1.e-9 ) ! Salt content drift (psu*km3)
355	CALL iom_put( 'bgvolssh' , zdiff_v1(tile_b) * 1.e-9 ) ! volume ssh drift (km3)
356	CALL iom_put( 'bgvole3t' , zdiff_v2(tile_b) * 1.e-9 ) ! volume e3t drift (km3)
357	!
358	! IF( lwp ) THEN
359	! WRITE(numout,*)
360	! WRITE(numout,*) 'dia_hsb : last time step hsb diagnostics: at it= ', kt,' date= ', ndastp
361	! WRITE(numout,*) '~~~~~~~'
362	! WRITE(numout,*) ' Temperature drift = ', zdiff_hc(tile_b) / zvol_tot(tile_b), ' C'
363	! WRITE(numout,*) ' Salinity12 drift = ', zdiff_sc(tile_b) / zvol_tot(tile_b), ' PSU'
364	! WRITE(numout,) ' volume ssh drift = ', zdiff_v1(tile_b) 1.e-9 , ' km^3'
365	! WRITE(numout,) ' volume e3t drift = ', zdiff_v2(tile_b) 1.e-9 , ' km^3'
366	! ENDIF
367	ELSE
368	CALL iom_put( 'bgtemper' , zdiff_hc1(tile_b) / zvol_tot(tile_b)) ! Heat content drift (C)
369	CALL iom_put( 'bgsaline' , zdiff_sc1(tile_b) / zvol_tot(tile_b)) ! Salt content drift (PSU)
370	CALL iom_put( 'bgheatco' , zdiff_hc1(tile_b) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
371	CALL iom_put( 'bgheatfx' , zdiff_hc1(tile_b) * rho0 * rcp / & ! Heat flux drift (W/m2)
372	& ( surf_tot * kts * rn_Dt ) )
373	CALL iom_put( 'bgsaltco' , zdiff_sc1(tile_b) * 1.e-9 ) ! Salt content drift (psu*km3)
374	CALL iom_put( 'bgvolssh' , zdiff_v1(tile_b) * 1.e-9 ) ! volume ssh drift (km3)
375	CALL iom_put( 'bgmistem' , zerr_hc1(tile_b) / zvol_tot(tile_b) ) ! hc - error due to free surface (C)
376	CALL iom_put( 'bgmissal' , zerr_sc1(tile_b) / zvol_tot(tile_b) ) ! sc - error due to free surface (psu)
377	ENDIF
378	!
379	IF( lrst_oce ) CALL dia_hsb_rst( kts, Kmm, tile_n, 'WRITE' )
380	!
381	END IF
382	IF ( kt == nitend ) THEN
383
384	zdiff_v1(tile_n) = zdiff_v1(tile_n) - frc_v(tile_n)
385	IF( .NOT.ln_linssh ) zdiff_v2(tile_n) = zdiff_v2(tile_n) - frc_v(tile_n)
386	zdiff_hc(tile_n) = zdiff_hc(tile_n) - frc_t(tile_n)
387	zdiff_sc(tile_n) = zdiff_sc(tile_n) - frc_s(tile_n)
388	IF( ln_linssh ) THEN
389	zdiff_hc1(tile_n) = zdiff_hc (tile_n) + z_ssh_hc(tile_n)
390	zdiff_sc1(tile_n) = zdiff_sc (tile_n) + z_ssh_sc(tile_n)
391	zerr_hc1 (tile_n) = z_ssh_hc(tile_n) - frc_wn_t(tile_n)
392	zerr_sc1 (tile_n) = z_ssh_sc(tile_n) - frc_wn_s(tile_n)
393	ENDIF
394	!!gm to be added ?
395	! IF( ln_linssh ) THEN ! fixed volume, add the ssh contribution
396	! zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * sshn(:,:) )
397	! ENDIF
398	!!gm end
399
400	CALL iom_put( 'bgfrcvol' , frc_v(tile_n) * 1.e-9 ) ! vol - surface forcing (km3)
401	CALL iom_put( 'bgfrctem' , frc_t(tile_n) * rho0 * rcp * 1.e-20 ) ! hc - surface forcing (1.e20 J)
402	CALL iom_put( 'bgfrchfx' , frc_t(tile_n) * rho0 * rcp / & ! hc - surface forcing (W/m2)
403	& ( surf_tot * kt * rn_Dt ) )
404	CALL iom_put( 'bgfrcsal' , frc_s(tile_n) * 1.e-9 ) ! sc - surface forcing (psu*km3)
405	IF( .NOT. ln_linssh ) THEN
406	CALL iom_put( 'bgtemper' , zdiff_hc(tile_n) / zvol_tot(tile_n) ) ! Temperature drift (C)
407	CALL iom_put( 'bgsaline' , zdiff_sc(tile_n) / zvol_tot(tile_n) ) ! Salinity drift (PSU)
408	CALL iom_put( 'bgheatco' , zdiff_hc(tile_n) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
409	CALL iom_put( 'bgheatfx' , zdiff_hc(tile_n) * rho0 * rcp / & ! Heat flux drift (W/m2)
410	& ( surf_tot * kt * rn_Dt ) )
411	CALL iom_put( 'bgsaltco' , zdiff_sc(tile_n) * 1.e-9 ) ! Salt content drift (psu*km3)
412	CALL iom_put( 'bgvolssh' , zdiff_v1(tile_n) * 1.e-9 ) ! volume ssh drift (km3)
413	CALL iom_put( 'bgvole3t' , zdiff_v2(tile_n) * 1.e-9 ) ! volume e3t drift (km3)
414	!
415	IF( kt == nitend .AND. lwp ) THEN
416	WRITE(numout,*)
417	WRITE(numout,*) 'dia_hsb : last time step hsb diagnostics: at it= ', kt,' date= ', ndastp
418	WRITE(numout,*) '~~~~~~~'
419	WRITE(numout,*) ' Temperature drift = ', zdiff_hc(tile_n) / zvol_tot(tile_n), ' C'
420	WRITE(numout,*) ' Salinity drift = ', zdiff_sc(tile_n) / zvol_tot(tile_n), ' PSU'
421	WRITE(numout,) ' volume ssh drift = ', zdiff_v1(tile_n) 1.e-9 , ' km^3'
422	WRITE(numout,) ' volume e3t drift = ', zdiff_v2(tile_n) 1.e-9 , ' km^3'
423	!
424	ENDIF
425	!
426	ELSE
427	CALL iom_put( 'bgtemper' , zdiff_hc1(tile_n) / zvol_tot(tile_n)) ! Heat content drift (C)
428	CALL iom_put( 'bgsaline' , zdiff_sc1(tile_n) / zvol_tot(tile_n)) ! Salt content drift (PSU)
429	CALL iom_put( 'bgheatco' , zdiff_hc1(tile_n) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
430	CALL iom_put( 'bgheatfx' , zdiff_hc1(tile_n) * rho0 * rcp / & ! Heat flux drift (W/m2)
431	& ( surf_tot * kt * rn_Dt ) )
432	CALL iom_put( 'bgsaltco' , zdiff_sc1(tile_n) * 1.e-9 ) ! Salt content drift (psu*km3)
433	CALL iom_put( 'bgvolssh' , zdiff_v1(tile_n) * 1.e-9 ) ! volume ssh drift (km3)
434	CALL iom_put( 'bgmistem' , zerr_hc1(tile_n) / zvol_tot(tile_n) ) ! hc - error due to free surface (C)
435	CALL iom_put( 'bgmissal' , zerr_sc1(tile_n) / zvol_tot(tile_n) ) ! sc - error due to free surface (psu)
436	ENDIF
437	!
438	!Last step, don't need restart
439	!IF( lrst_oce ) CALL dia_hsb_rst( kts, Kmm, tile_n, 'WRITE' )
440	!
441	END IF
442	IF( ln_timing ) CALL timing_stop('dia_hsb')
443	!
444	END SUBROUTINE dia_hsb
445
446
447	SUBROUTINE dia_hsb_rst( kt, Kmm, tile, cdrw )
448	!!---------------------------------------------------------------------
449	!! * ROUTINE dia_hsb_rst *
450	!!
451	!! ** Purpose : Read or write DIA file in restart file
452	!!
453	!! ** Method : use of IOM library
454	!!----------------------------------------------------------------------
455	INTEGER , INTENT(in) :: kt ! ocean time-step
456	INTEGER , INTENT(in) :: Kmm ! ocean time level index
457	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
458	INTEGER , INTENT(in) :: tile ! host tile
459	!
460	INTEGER :: ji, jj, jk ! dummy loop indices
461	!!----------------------------------------------------------------------
462	!
463	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
464	IF( ln_rstart ) THEN !* Read the restart file
465	!
466	IF(lwp) WRITE(numout,*)
467	IF(lwp) WRITE(numout,*) ' dia_hsb_rst : read hsb restart at it= ', kt,' date= ', ndastp
468	IF(lwp) WRITE(numout,*)
469	CALL iom_get( numror, 'frc_v', frc_v(tile) )
470	CALL iom_get( numror, 'frc_t', frc_t(tile) )
471	CALL iom_get( numror, 'frc_s', frc_s(tile) )
472	IF( ln_linssh ) THEN
473	CALL iom_get( numror, 'frc_wn_t', frc_wn_t(tile) )
474	CALL iom_get( numror, 'frc_wn_s', frc_wn_s(tile) )
475	ENDIF
476	CALL iom_get( numror, jpdom_auto, 'surf_ini' , surf_ini ) ! ice sheet coupling
477	CALL iom_get( numror, jpdom_auto, 'ssh_ini' , ssh_ini )
478	CALL iom_get( numror, jpdom_auto, 'e3t_ini' , e3t_ini )
479	CALL iom_get( numror, jpdom_auto, 'tmask_ini' , tmask_ini )
480	CALL iom_get( numror, jpdom_auto, 'hc_loc_ini', hc_loc_ini )
481	CALL iom_get( numror, jpdom_auto, 'sc_loc_ini', sc_loc_ini )
482	IF( ln_linssh ) THEN
483	CALL iom_get( numror, jpdom_auto, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
484	CALL iom_get( numror, jpdom_auto, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
485	ENDIF
486	ELSE
487	IF(lwp) WRITE(numout,*)
488	IF(lwp) WRITE(numout,*) ' dia_hsb_rst : initialise hsb at initial state '
489	IF(lwp) WRITE(numout,*)
490	surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:) ! initial ocean surface
491	ssh_ini(:,:) = ssh(:,:,Kmm) ! initial ssh
492	DO jk = 1, jpk
493	! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
494	e3t_ini (:,:,jk) = e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial vertical scale factors
495	tmask_ini (:,:,jk) = tmask(:,:,jk) ! initial mask
496	hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial heat content
497	sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial salt content
498	END DO
499	d_surf_ini = surf_ini
500	d_e3t_ini = e3t_ini
501	d_tmask_ini = tmask_ini
502	d_hc_loc_ini = hc_loc_ini
503	d_sc_loc_ini = sc_loc_ini
504	frc_v(tile) = 0._wp ! volume trend due to forcing
505	frc_t(tile) = 0._wp ! heat content - - - -
506	frc_s(tile) = 0._wp ! salt content - - - -
507	IF( ln_linssh ) THEN
508	IF( ln_isfcav ) THEN
509	DO ji = 1, jpi
510	DO jj = 1, jpj
511	ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) ! initial heat content in ssh
512	ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) ! initial salt content in ssh
513	END DO
514	END DO
515	ELSE
516	ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) ! initial heat content in ssh
517	ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) ! initial salt content in ssh
518	END IF
519	frc_wn_t(tile) = 0._wp ! initial heat content misfit due to free surface
520	frc_wn_s(tile) = 0._wp ! initial salt content misfit due to free surface
521	ENDIF
522	ENDIF
523	!
524	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
525	! ! -------------------
526	IF(lwp) WRITE(numout,*)
527	IF(lwp) WRITE(numout,*) ' dia_hsb_rst : write restart at it= ', kt,' date= ', ndastp
528	IF(lwp) WRITE(numout,*)
529	!
530	CALL iom_rstput( kt, nitrst, numrow, 'frc_v', frc_v(tile) )
531	CALL iom_rstput( kt, nitrst, numrow, 'frc_t', frc_t(tile) )
532	CALL iom_rstput( kt, nitrst, numrow, 'frc_s', frc_s(tile) )
533	IF( ln_linssh ) THEN
534	CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t(tile) )
535	CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s(tile) )
536	ENDIF
537	CALL iom_rstput( kt, nitrst, numrow, 'surf_ini' , surf_ini ) ! ice sheet coupling
538	CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini' , ssh_ini )
539	CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini' , e3t_ini )
540	CALL iom_rstput( kt, nitrst, numrow, 'tmask_ini' , tmask_ini )
541	CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini )
542	CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
543	IF( ln_linssh ) THEN
544	CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
545	CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
546	ENDIF
547	!
548	ENDIF
549	!
550	END SUBROUTINE dia_hsb_rst
551
552
553	SUBROUTINE dia_hsb_init( Kmm )
554	!!---------------------------------------------------------------------------
555	!! * ROUTINE dia_hsb *
556	!!
557	!! ** Purpose: Initialization for the heat salt volume budgets
558	!!
559	!! ** Method : Compute initial heat content, salt content and volume
560	!!
561	!! ** Action : - Compute initial heat content, salt content and volume
562	!! - Initialize forcing trends
563	!! - Compute coefficients for conversion
564	!!---------------------------------------------------------------------------
565	INTEGER, INTENT(in) :: Kmm ! time level index
566	!
567	INTEGER :: ierror, ios ! local integer
568	INTEGER :: i, istat ! local integer
569	!!
570	NAMELIST/namhsb/ ln_diahsb
571	!!----------------------------------------------------------------------
572	!
573	IF(lwp) THEN
574	WRITE(numout,*)
575	WRITE(numout,*) 'dia_hsb_init : heat and salt budgets diagnostics'
576	WRITE(numout,*) '~~~~~~~~~~~~ '
577	ENDIF
578	READ ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
579	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist' )
580	READ ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
581	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist' )
582	IF(lwm) WRITE( numond, namhsb )
583
584	IF(lwp) THEN
585	WRITE(numout,*) ' Namelist namhsb :'
586	WRITE(numout,*) ' check the heat and salt budgets (T) or not (F) ln_diahsb = ', ln_diahsb
587	ENDIF
588	!
589	IF( .NOT. ln_diahsb ) RETURN
590
591	! ------------------- !
592	! 1 - Allocate memory !
593	! ------------------- !
594
595	CALL setdevice()
596	!Device data associate to PUBLIC arrays
597	ALLOCATE(d_e3t (jpi,jpj,jpk,jpt) ) !
598	ALLOCATE(d_tmask (jpi,jpj,jpk) ) !
599	ALLOCATE(d_tmask_ini (jpi,jpj,jpk) ) !
600	ALLOCATE(d_tmask_h (jpi,jpj) ) !
601	ALLOCATE(d_ts (jpi,jpj,jpk,2,jpj) ) !
602	!Device data associate to LOCAL/DEVICE arrays !
603	ALLOCATE(d_surf (jpi,jpj) ) !
604	ALLOCATE(d_surf_ini (jpi,jpj) ) !
605	ALLOCATE(d_hc_loc_ini (jpi,jpj,jpk) ) !
606	ALLOCATE(d_sc_loc_ini (jpi,jpj,jpk) ) !
607	ALLOCATE(d_e3t_ini (jpi,jpj,jpk) ) !
608	ALLOCATE(d_zwrkv (jpi,jpj,jpkm1) ) !
609	ALLOCATE(d_zwrkh (jpi,jpj,jpkm1) ) !
610	ALLOCATE(d_zwrks (jpi,jpj,jpkm1) ) !
611	ALLOCATE(d_zwrk (jpi,jpj,jpkm1) ) !
612	ALLOCATE(h_ztmpv(2),h_ztmph(2),h_ztmps(2),h_ztmp(2)) !
613
614	DO i = 1, nstreams !Create Streams
615	istat = cudaStreamCreate(stream(i))
616	IF( istat /= 0 ) THEN
617	CALL ctl_stop( 'dia_hsb_init: error in Stream creation' ) ; RETURN
618	ENDIF
619	END DO
620	!
621	!Pinned reallocation step non constant
622	istat = cudaHostRegister(C_LOC(ts ), sizeof(ts ), cudaHostRegisterMapped)
623	istat = cudaHostRegister(C_LOC(e3t), sizeof(e3t), cudaHostRegisterMapped)
624	IF( istat /= 0 ) THEN
625	CALL ctl_stop( 'dia_hsb_init: unable to pin host memory to GPU' ) ; RETURN
626	ENDIF
627
628	ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
629	& e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), tmask_ini(jpi,jpj,jpk),STAT=ierror )
630	IF( ierror > 0 ) THEN
631	CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' ) ; RETURN
632	ENDIF
633
634	IF( ln_linssh ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
635	IF( ierror > 0 ) THEN
636	CALL ctl_stop( 'dia_hsb: unable to allocate ssh_hc_loc_ini' ) ; RETURN
637	ENDIF
638
639	! ----------------------------------------------- !
640	! 2 - Time independant variables and file opening !
641	! ----------------------------------------------- !
642	surf(:,:) = e1e2t(:,:) * tmask_i(:,:) ! masked surface grid cell area
643	surf_tot = glob_sum( 'diahsb', surf(:,:) ) ! total ocean surface area
644
645	d_surf = surf
646	d_surf_ini = surf_ini
647	d_e3t_ini = e3t_ini
648	d_tmask = tmask
649	d_tmask_ini = tmask_ini
650	d_tmask_h = tmask_h
651	h_ztmp = 0.0
652
653	IF( ln_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )
654	!
655	! ---------------------------------- !
656	! 4 - initial conservation variables !
657	! ---------------------------------- !
658	CALL dia_hsb_rst( nit000, Kmm, 1, 'READ' ) !* read or initialize all required files
659
660
661
662	!
663	END SUBROUTINE dia_hsb_init
664
665	!!======================================================================
666	END MODULE diahsb

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format

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

Context Navigation

source: NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/diahsb_gpu.h90 @ 15063

Download in other formats: