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.

diaar5.F90 in NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA/diaar5.F90 @ 12532

Last change on this file since 12532 was 12532, checked in by davestorkey, 4 years ago
UKMO/NEMO_4.0.1_biharmonic_GM: science changes.
File size: 16.9 KB

Line
1	MODULE diaar5
2	!!======================================================================
3	!! * MODULE diaar5 *
4	!! AR5 diagnostics
5	!!======================================================================
6	!! History : 3.2 ! 2009-11 (S. Masson) Original code
7	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
8	!!----------------------------------------------------------------------
9	!! dia_ar5 : AR5 diagnostics
10	!! dia_ar5_init : initialisation of AR5 diagnostics
11	!!----------------------------------------------------------------------
12	USE oce ! ocean dynamics and active tracers
13	USE dom_oce ! ocean space and time domain
14	USE eosbn2 ! equation of state (eos_bn2 routine)
15	USE phycst ! physical constant
16	USE in_out_manager ! I/O manager
17	USE zdfddm
18	USE zdf_oce
19	!
20	USE lib_mpp ! distribued memory computing library
21	USE iom ! I/O manager library
22	USE fldread ! type FLD_N
23	USE timing ! preformance summary
24
25	IMPLICIT NONE
26	PRIVATE
27
28	PUBLIC dia_ar5 ! routine called in step.F90 module
29	PUBLIC dia_ar5_alloc ! routine called in nemogcm.F90 module
30	PUBLIC dia_ar5_hst ! heat/salt transport
31
32	REAL(wp) :: vol0 ! ocean volume (interior domain)
33	REAL(wp) :: area_tot ! total ocean surface (interior domain)
34	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: area ! cell surface (interior domain)
35	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: thick0 ! ocean thickness (interior domain)
36	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sn0 ! initial salinity
37
38	LOGICAL :: l_ar5
39
40	!! * Substitutions
41	# include "vectopt_loop_substitute.h90"
42	!!----------------------------------------------------------------------
43	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
44	!! $Id: diaar5.F90 10425 2018-12-19 21:54:16Z smasson $
45	!! Software governed by the CeCILL license (see ./LICENSE)
46	!!----------------------------------------------------------------------
47	CONTAINS
48
49	FUNCTION dia_ar5_alloc()
50	!!----------------------------------------------------------------------
51	!! * ROUTINE dia_ar5_alloc *
52	!!----------------------------------------------------------------------
53	INTEGER :: dia_ar5_alloc
54	!!----------------------------------------------------------------------
55	!
56	ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc )
57	!
58	CALL mpp_sum ( 'diaar5', dia_ar5_alloc )
59	IF( dia_ar5_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dia_ar5_alloc: failed to allocate arrays' )
60	!
61	END FUNCTION dia_ar5_alloc
62
63
64	SUBROUTINE dia_ar5( kt )
65	!!----------------------------------------------------------------------
66	!! * ROUTINE dia_ar5 *
67	!!
68	!! ** Purpose : compute and output some AR5 diagnostics
69	!!----------------------------------------------------------------------
70	!
71	INTEGER, INTENT( in ) :: kt ! ocean time-step index
72	!
73	INTEGER :: ji, jj, jk ! dummy loop arguments
74	REAL(wp) :: zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass
75	REAL(wp) :: zaw, zbw, zrw
76	!
77	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zarea_ssh , zbotpres ! 2D workspace
78	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpe ! 2D workspace
79	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zrhd , zrhop ! 3D workspace
80	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn ! 4D workspace
81
82	!!--------------------------------------------------------------------
83	IF( ln_timing ) CALL timing_start('dia_ar5')
84
85	IF( kt == nit000 ) CALL dia_ar5_init
86
87	IF( l_ar5 ) THEN
88	ALLOCATE( zarea_ssh(jpi,jpj) , zbotpres(jpi,jpj) )
89	ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) )
90	ALLOCATE( ztsn(jpi,jpj,jpk,jpts) )
91	zarea_ssh(:,:) = area(:,:) * sshn(:,:)
92	ENDIF
93	!
94	IF( iom_use( 'voltot' ) .OR. iom_use( 'sshtot' ) .OR. iom_use( 'sshdyn' ) ) THEN
95	! ! total volume of liquid seawater
96	zvolssh = SUM( zarea_ssh(:,:) )
97	CALL mpp_sum( 'diaar5', zvolssh )
98	zvol = vol0 + zvolssh
99
100	CALL iom_put( 'voltot', zvol )
101	CALL iom_put( 'sshtot', zvolssh / area_tot )
102	CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
103	!
104	ENDIF
105
106	IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' ) .OR. iom_use( 'sshsteric' ) ) THEN
107	!
108	ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) ! thermosteric ssh
109	ztsn(:,:,:,jp_sal) = sn0(:,:,:)
110	CALL eos( ztsn, zrhd, gdept_n(:,:,:) ) ! now in situ density using initial salinity
111	!
112	zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice
113	DO jk = 1, jpkm1
114	zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
115	END DO
116	IF( ln_linssh ) THEN
117	IF( ln_isfcav ) THEN
118	DO ji = 1, jpi
119	DO jj = 1, jpj
120	zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
121	END DO
122	END DO
123	ELSE
124	zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
125	END IF
126	!!gm
127	!!gm riceload should be added in both ln_linssh=T or F, no?
128	!!gm
129	END IF
130	!
131	zarho = SUM( area(:,:) * zbotpres(:,:) )
132	CALL mpp_sum( 'diaar5', zarho )
133	zssh_steric = - zarho / area_tot
134	CALL iom_put( 'sshthster', zssh_steric )
135
136	! ! steric sea surface height
137	CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) ) ! now in situ and potential density
138	!CW: commenting out the following line, as seems weird to zero the array before output & only way to get rhop output
139	!zrhop(:,:,jpk) = 0._wp
140	CALL iom_put( 'rhop', zrhop )
141	!
142	zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice
143	DO jk = 1, jpkm1
144	zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
145	END DO
146	IF( ln_linssh ) THEN
147	IF ( ln_isfcav ) THEN
148	DO ji = 1,jpi
149	DO jj = 1,jpj
150	zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
151	END DO
152	END DO
153	ELSE
154	zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
155	END IF
156	END IF
157	!
158	zarho = SUM( area(:,:) * zbotpres(:,:) )
159	CALL mpp_sum( 'diaar5', zarho )
160	zssh_steric = - zarho / area_tot
161	CALL iom_put( 'sshsteric', zssh_steric )
162
163	! ! ocean bottom pressure
164	zztmp = rau0 * grav * 1.e-4_wp ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
165	zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
166	CALL iom_put( 'botpres', zbotpres )
167	!
168	ENDIF
169
170	IF( iom_use( 'masstot' ) .OR. iom_use( 'temptot' ) .OR. iom_use( 'saltot' ) ) THEN
171	! ! Mean density anomalie, temperature and salinity
172	ztemp = 0._wp
173	zsal = 0._wp
174	DO jk = 1, jpkm1
175	DO jj = 1, jpj
176	DO ji = 1, jpi
177	zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
178	ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
179	zsal = zsal + zztmp * tsn(ji,jj,jk,jp_sal)
180	END DO
181	END DO
182	END DO
183	IF( ln_linssh ) THEN
184	IF( ln_isfcav ) THEN
185	DO ji = 1, jpi
186	DO jj = 1, jpj
187	ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)
188	zsal = zsal + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)
189	END DO
190	END DO
191	ELSE
192	ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
193	zsal = zsal + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
194	END IF
195	ENDIF
196	IF( lk_mpp ) THEN
197	CALL mpp_sum( 'diaar5', ztemp )
198	CALL mpp_sum( 'diaar5', zsal )
199	END IF
200	!
201	zmass = rau0 * ( zarho + zvol ) ! total mass of liquid seawater
202	ztemp = ztemp / zvol ! potential temperature in liquid seawater
203	zsal = zsal / zvol ! Salinity of liquid seawater
204	!
205	CALL iom_put( 'masstot', zmass )
206	CALL iom_put( 'temptot', ztemp )
207	CALL iom_put( 'saltot' , zsal )
208	!
209	ENDIF
210
211	IF( iom_use( 'tnpeo' )) THEN
212	! Work done against stratification by vertical mixing
213	! Exclude points where rn2 is negative as convection kicks in here and
214	! work is not being done against stratification
215	ALLOCATE( zpe(jpi,jpj) )
216	zpe(:,:) = 0._wp
217	IF( ln_zdfddm ) THEN
218	DO jk = 2, jpk
219	DO jj = 1, jpj
220	DO ji = 1, jpi
221	IF( rn2(ji,jj,jk) > 0._wp ) THEN
222	zrw = ( gdepw_n(ji,jj,jk ) - gdept_n(ji,jj,jk) ) &
223	& / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
224	!!gm this can be reduced to : (depw-dept) / e3w (NB idem dans bn2 !)
225	! zrw = ( gdept_n(ji,jj,jk) - gdepw_n(ji,jj,jk) ) / e3w_n(ji,jj,jk)
226	!!gm end
227	!
228	zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
229	zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
230	!
231	zpe(ji, jj) = zpe(ji, jj) &
232	& - grav * ( avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) &
233	& - avs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
234	ENDIF
235	END DO
236	END DO
237	END DO
238	ELSE
239	DO jk = 1, jpk
240	DO ji = 1, jpi
241	DO jj = 1, jpj
242	zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
243	END DO
244	END DO
245	END DO
246	ENDIF
247	!!gm useless lbc_lnk since the computation above is performed over 1:jpi & 1:jpj
248	!!gm CALL lbc_lnk( 'diaar5', zpe, 'T', 1._wp)
249	CALL iom_put( 'tnpeo', zpe )
250	DEALLOCATE( zpe )
251	ENDIF
252
253	IF( l_ar5 ) THEN
254	DEALLOCATE( zarea_ssh , zbotpres )
255	DEALLOCATE( zrhd , zrhop )
256	DEALLOCATE( ztsn )
257	ENDIF
258	!
259	IF( ln_timing ) CALL timing_stop('dia_ar5')
260	!
261	END SUBROUTINE dia_ar5
262
263
264	SUBROUTINE dia_ar5_hst( ktra, cptr, pua, pva )
265	!!----------------------------------------------------------------------
266	!! * ROUTINE dia_ar5_htr *
267	!!----------------------------------------------------------------------
268	!! Wrapper for heat transport calculations
269	!! Called from all advection and/or diffusion routines
270	!!----------------------------------------------------------------------
271	INTEGER , INTENT(in ) :: ktra ! tracer index
272	CHARACTER(len=3) , INTENT(in) :: cptr ! transport type 'adv'/'ldf'
273	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pua ! 3D input array of advection/diffusion
274	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pva ! 3D input array of advection/diffusion
275	!
276	INTEGER :: ji, jj, jk
277	REAL(wp), DIMENSION(jpi,jpj) :: z2d
278
279
280	z2d(:,:) = pua(:,:,1)
281	DO jk = 1, jpkm1
282	DO jj = 2, jpjm1
283	DO ji = fs_2, fs_jpim1 ! vector opt.
284	z2d(ji,jj) = z2d(ji,jj) + pua(ji,jj,jk)
285	END DO
286	END DO
287	END DO
288	CALL lbc_lnk( 'diaar5', z2d, 'U', -1. )
289	IF( cptr == 'adv' ) THEN
290	IF( ktra == jp_tem ) CALL iom_put( "uadv_heattr" , rau0_rcp * z2d ) ! advective heat transport in i-direction
291	IF( ktra == jp_sal ) CALL iom_put( "uadv_salttr" , rau0 * z2d ) ! advective salt transport in i-direction
292	ENDIF
293	IF( cptr == 'ldf' ) THEN
294	IF( ktra == jp_tem ) CALL iom_put( "udiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in i-direction
295	IF( ktra == jp_sal ) CALL iom_put( "udiff_salttr" , rau0 * z2d ) ! diffusive salt transport in i-direction
296	ENDIF
297	!
298	z2d(:,:) = pva(:,:,1)
299	DO jk = 1, jpkm1
300	DO jj = 2, jpjm1
301	DO ji = fs_2, fs_jpim1 ! vector opt.
302	z2d(ji,jj) = z2d(ji,jj) + pva(ji,jj,jk)
303	END DO
304	END DO
305	END DO
306	CALL lbc_lnk( 'diaar5', z2d, 'V', -1. )
307	IF( cptr == 'adv' ) THEN
308	IF( ktra == jp_tem ) CALL iom_put( "vadv_heattr" , rau0_rcp * z2d ) ! advective heat transport in j-direction
309	IF( ktra == jp_sal ) CALL iom_put( "vadv_salttr" , rau0 * z2d ) ! advective salt transport in j-direction
310	ENDIF
311	IF( cptr == 'ldf' ) THEN
312	IF( ktra == jp_tem ) CALL iom_put( "vdiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in j-direction
313	IF( ktra == jp_sal ) CALL iom_put( "vdiff_salttr" , rau0 * z2d ) ! diffusive salt transport in j-direction
314	ENDIF
315
316	END SUBROUTINE dia_ar5_hst
317
318
319	SUBROUTINE dia_ar5_init
320	!!----------------------------------------------------------------------
321	!! * ROUTINE dia_ar5_init *
322	!!
323	!! ** Purpose : initialization for AR5 diagnostic computation
324	!!----------------------------------------------------------------------
325	INTEGER :: inum
326	INTEGER :: ik
327	INTEGER :: ji, jj, jk ! dummy loop indices
328	REAL(wp) :: zztmp
329	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: zsaldta ! Jan/Dec levitus salinity
330	!
331	!!----------------------------------------------------------------------
332	!
333	l_ar5 = .FALSE.
334	IF( iom_use( 'voltot' ) .OR. iom_use( 'sshtot' ) .OR. iom_use( 'sshdyn' ) .OR. &
335	& iom_use( 'masstot' ) .OR. iom_use( 'temptot' ) .OR. iom_use( 'saltot' ) .OR. &
336	& iom_use( 'botpres' ) .OR. iom_use( 'sshthster' ) .OR. iom_use( 'sshsteric' ) ) L_ar5 = .TRUE.
337
338	IF( l_ar5 ) THEN
339	!
340	! ! allocate dia_ar5 arrays
341	IF( dia_ar5_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
342
343	area(:,:) = e1e2t(:,:) * tmask_i(:,:)
344
345	area_tot = SUM( area(:,:) ) ; CALL mpp_sum( 'diaar5', area_tot )
346
347	vol0 = 0._wp
348	thick0(:,:) = 0._wp
349	DO jk = 1, jpkm1
350	vol0 = vol0 + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
351	thick0(:,:) = thick0(:,:) + tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
352	END DO
353	CALL mpp_sum( 'diaar5', vol0 )
354
355	IF( iom_use( 'sshthster' ) ) THEN
356	ALLOCATE( zsaldta(jpi,jpj,jpj,jpts) )
357	CALL iom_open ( 'sali_ref_clim_monthly', inum )
358	CALL iom_get ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1 )
359	CALL iom_get ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 )
360	CALL iom_close( inum )
361
362	sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )
363	sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
364	IF( ln_zps ) THEN ! z-coord. partial steps
365	DO jj = 1, jpj ! interpolation of salinity at the last ocean level (i.e. the partial step)
366	DO ji = 1, jpi
367	ik = mbkt(ji,jj)
368	IF( ik > 1 ) THEN
369	zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
370	sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
371	ENDIF
372	END DO
373	END DO
374	ENDIF
375	!
376	DEALLOCATE( zsaldta )
377	ENDIF
378	!
379	ENDIF
380	!
381	END SUBROUTINE dia_ar5_init
382
383	!!======================================================================
384	END MODULE diaar5

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

Download in other formats: