Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

diaregmean.F90

Last change on this file was 11085, checked in by rrenshaw, 5 years ago
code fixes
File size: 52.9 KB

Line
1	MODULE diaregmean
2	!!======================================================================
3	!! * MODULE diaharm *
4	!! Timeseries of Regional Means
5	!!======================================================================
6	!! History : 3.6 ! 11/2016 (J Tinker) Original code
7	!!----------------------------------------------------------------------
8	USE oce ! ocean dynamics and tracers variables
9	USE dom_oce ! ocean space and time domain
10	USE in_out_manager ! I/O units
11	USE iom ! I/0 library
12	USE wrk_nemo ! working arrays
13	USE diapea ! PEA
14	USE zdfmxl ! MLD
15	USE sbc_oce
16	#if defined key_diaar5
17	USE diaar5
18	#endif
19
20
21	#if defined key_fabm
22	USE trc
23	USE par_fabm
24	#endif
25
26	IMPLICIT NONE
27	PRIVATE
28
29	LOGICAL , PUBLIC :: ln_diaregmean ! region mean calculation
30	PUBLIC dia_regmean_init ! routine called by nemogcm.F90
31	PUBLIC dia_regmean ! routine called by diawri.F90
32	PUBLIC dia_calctmb_region_mean ! routine called by diatmb.F90
33
34
35
36	LOGICAL :: ln_diaregmean_ascii ! region mean calculation ascii output
37	LOGICAL :: ln_diaregmean_bin ! region mean calculation binary output
38	LOGICAL :: ln_diaregmean_nc ! region mean calculation netcdf output
39	LOGICAL :: ln_diaregmean_diaar5 ! region mean calculation including AR5 SLR terms
40	LOGICAL :: ln_diaregmean_diasbc ! region mean calculation including Surface BC
41	LOGICAL :: ln_diaregmean_karamld ! region mean calculation including kara mld terms
42	LOGICAL :: ln_diaregmean_pea ! region mean calculation including pea terms
43
44
45	LOGICAL :: ln_diaregmean_bgc ! region mean calculation including BGC terms
46
47
48
49	REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: tmp_region_mask_real ! tempory region_mask of reals
50	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: region_mask ! region_mask matrix
51	INTEGER :: nmasks ! Number of mask files in region_mask.nc file -
52	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: nreg_mat ! Number of regions in each mask
53
54	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_field_mat !: temporary region_mask
55	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_field_AR5_mat !: temporary region_mask
56	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_field_SBC_mat !: temporary region_mask
57	INTEGER :: tmp_field_cnt ! tmp_field_cnt integer
58	!!----------------------------------------------------------------------
59	!! NEMO/OPA 3.6 , NEMO Consortium (2014)
60	!! $Id$
61	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
62	!!----------------------------------------------------------------------
63	CONTAINS
64
65	SUBROUTINE dia_regmean_init
66	!!---------------------------------------------------------------------------
67	!! * ROUTINE dia_regmean_init *
68	!!
69	!! ** Purpose: Initialization of region mask namelist
70	!!
71	!! ** Method : Read namelist
72	!! History
73	!! 3.6 ! 11-16 (J Tinker) Routine to initialize dia_regmean
74	!!---------------------------------------------------------------------------
75	!!
76	INTEGER :: ios ! Local integer output status for namelist read
77	INTEGER :: inum ! temporary logical unit ! copied from DOM/domzgr.F90
78	INTEGER :: ierr ! error integer for IOM_get
79	INTEGER :: idmaskvar ! output of iom_varid
80	INTEGER :: maskno ! counter for number of masks
81	INTEGER :: jj,ji ! i and j index
82	INTEGER :: tmpint ! temporary integer
83	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tmpregion !: temporary region_mask
84	INTEGER, DIMENSION(3) :: zdimsz ! number of elements in each of the 3 dimensions (i.e., lon, lat, no of masks, 297, 375, 4) for an array
85	INTEGER :: zndims ! number of dimensions in an array (i.e. 3, )
86
87
88	#if defined key_fabm
89	INTEGER :: js,jl,jn, tmp_dummy
90
91	CHARACTER (len=120) :: tmp_name,tmp_long_name, tmp_unit
92
93	INTEGER :: BGC_nlevs,nBGC_output, bgci
94	CHARACTER(len = 10), ALLOCATABLE, DIMENSION(:) :: BGC_stat_name(:),BGC_lev_name(:),BGC_output_var(:)
95	#endif
96
97	!
98	NAMELIST/nam_diaregmean/ ln_diaregmean,ln_diaregmean_ascii,ln_diaregmean_bin,ln_diaregmean_nc,&
99	& ln_diaregmean_karamld, ln_diaregmean_pea,ln_diaregmean_diaar5,ln_diaregmean_diasbc,ln_diaregmean_bgc
100
101
102	! read in Namelist.
103	!!----------------------------------------------------------------------
104	!
105	! Defaults are generally negative, regional means output switched off
106	ln_diaregmean = .FALSE. ! don't write regional means
107	ln_diaregmean_bgc = .FALSE. ! don't write biogeochemistry regional means
108	ln_diaregmean_karamld = .FALSE. ! don't write mixed layer depth
109	ln_diaregmean_pea = .FALSE. ! don't write potential energy anomaly
110	ln_diaregmean_diaar5 = .FALSE. ! don't write AR5 SLR terms
111	ln_diaregmean_diasbc = .FALSE. ! don't write surface forcing terms
112
113	REWIND ( numnam_ref ) ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
114	READ ( numnam_ref, nam_diaregmean, IOSTAT=ios, ERR= 901 )
115	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in reference namelist', lwp )
116
117	REWIND( numnam_cfg ) ! Namelist nam_diatmb in configuration namelist TMB diagnostics
118	READ ( numnam_cfg, nam_diaregmean, IOSTAT = ios, ERR = 902 )
119	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in configuration namelist', lwp )
120	IF(lwm) WRITE ( numond, nam_diaregmean )
121
122	IF(lwp) THEN ! Control print
123	WRITE(numout,*)
124	WRITE(numout,*) 'dia_regmean_init : Output regional mean Diagnostics'
125	WRITE(numout,*) '~~~~~~~~~~~~'
126	WRITE(numout,*) 'Namelist nam_regmean : set regmeanoutputs '
127	WRITE(numout,*) 'Switch for regmean diagnostics (T) or not (F) ln_diaregmean = ', ln_diaregmean
128	WRITE(numout,*) 'Switch for regmean ascii output (T) or not (F) ln_diaregmean_ascii = ', ln_diaregmean_ascii
129	WRITE(numout,*) 'Switch for regmean binary output (T) or not (F) ln_diaregmean_bin = ', ln_diaregmean_bin
130	WRITE(numout,*) 'Switch for regmean netcdf output (T) or not (F) ln_diaregmean_nc = ', ln_diaregmean_nc
131	WRITE(numout,*) 'Switch for regmean kara mld terms (T) or not (F) ln_diaregmean_karamld = ', ln_diaregmean_karamld
132	WRITE(numout,*) 'Switch for regmean PEA terms (T) or not (F) ln_diaregmean_pea = ', ln_diaregmean_pea
133	WRITE(numout,*) 'Switch for regmean AR5 SLR terms (T) or not (F) ln_diaregmean_diaar5 = ', ln_diaregmean_diaar5
134	WRITE(numout,*) 'Switch for regmean Surface forcing terms (T) or not (F) ln_diaregmean_diasbc = ', ln_diaregmean_diasbc
135	WRITE(numout,*) 'Switch for regmean BioGeoChemistry terms (T) or not (F) ln_diaregmean_bgc = ', ln_diaregmean_bgc
136	ENDIF
137
138
139	ALLOCATE( tmp_field_mat(jpi,jpj,19), STAT= ierr ) !SS/NB/DT/ZA/VA T/S, SSH, MLD, PEA, PEAT, PEAS
140	IF( ierr /= 0 ) CALL ctl_stop( 'tmp_field_mat: failed to allocate tmp_field_mat array' )
141	tmp_field_mat(:,:,:) = 0.
142	tmp_field_cnt = 0
143
144	IF(ln_diaregmean_diaar5) THEN
145	ALLOCATE( tmp_field_AR5_mat(jpi,jpj,4), STAT= ierr ) !SLR terms
146	IF( ierr /= 0 ) CALL ctl_stop( 'tmp_field_AR5_mat: failed to allocate tmp_field_AR5_mat array' )
147	tmp_field_AR5_mat(:,:,:) = 0.
148	ENDIF
149
150	IF(ln_diaregmean_diasbc) THEN
151	ALLOCATE( tmp_field_SBC_mat(jpi,jpj,7), STAT= ierr ) !SBC terms
152	IF( ierr /= 0 ) CALL ctl_stop( 'tmp_field_SBC_mat: failed to allocate tmp_field_SBC_mat array' )
153	tmp_field_SBC_mat(:,:,:) = 0.
154	ENDIF
155
156
157	#if defined key_fabm
158	! as there are so many BGC variables, write out the necessary iodef.xml and field_def.xml entries into ocean.output
159
160	IF(ln_diaregmean_bgc) THEN
161	IF(lwp) THEN ! Control print
162
163	BGC_nlevs = 5
164	ALLOCATE( BGC_stat_name(6),BGC_lev_name(BGC_nlevs))
165	nBGC_output = 16
166	ALLOCATE( BGC_output_var(nBGC_output))
167
168	BGC_output_var(1) = 'N1_p'
169	BGC_output_var(2) = 'N3_n'
170	BGC_output_var(3) = 'N4_n'
171	BGC_output_var(4) = 'N5_s'
172	BGC_output_var(5) = 'O2_o'
173	BGC_output_var(6) = 'P1_Chl'
174	BGC_output_var(7) = 'P2_Chl'
175	BGC_output_var(8) = 'P3_Chl'
176	BGC_output_var(9) = 'P4_Chl'
177	BGC_output_var(10) = 'P1_c'
178	BGC_output_var(11) = 'P2_c'
179	BGC_output_var(12) = 'P3_c'
180	BGC_output_var(13) = 'P4_c'
181	BGC_output_var(14) = 'Z4_c'
182	BGC_output_var(15) = 'Z5_c'
183	BGC_output_var(16) = 'Z6_c'
184
185	BGC_stat_name(1) = '_ave'
186	BGC_stat_name(2) = '_tot'
187	BGC_stat_name(3) = '_var'
188	BGC_stat_name(4) = '_cnt'
189	BGC_stat_name(5) = '_reg_id'
190	BGC_stat_name(6) = '_mask_id'
191	BGC_lev_name(1) = 'top'
192	BGC_lev_name(2) = 'bot'
193	BGC_lev_name(3) = 'dif'
194	BGC_lev_name(4) = 'zav'
195	BGC_lev_name(5) = 'vol'
196
197
198	WRITE(numout,*) ''
199	WRITE(numout,*) 'diaregmean BGC field_def.xml entries'
200	WRITE(numout,*) ''
201
202
203	DO jn=1,jp_fabm ! State loop
204	DO js=1,6
205	DO jl=1,BGC_nlevs
206
207	tmp_name=TRIM( TRIM("reg_")//TRIM(BGC_lev_name(jl))//TRIM("_")//TRIM(ctrcnm(jn))// TRIM(BGC_stat_name(js)) )
208
209	tmp_long_name = TRIM(ctrcln(jn))
210	tmp_unit = TRIM(ctrcun(jn))
211
212	! Where using volume integrated values, change units...
213
214	IF ((jl .EQ. 5) .AND. (js .EQ. 2)) then
215	SELECT CASE (trim(tmp_unit))
216	CASE ('mg C/m^3') ; tmp_unit = 'Mg C (T C)' !'mg C/m^3'
217	CASE ('mg/m^3') ; tmp_unit = 'Mg (T)' !'mg/m^3'
218	CASE ('mmol C/m^3') ; tmp_unit = 'Mmol C' !'mmol C/m^3'
219	CASE ('mmol N/m^3') ; tmp_unit = 'Mmol N' !'mmol N/m^3'
220	CASE ('mmol O_2/m^3') ; tmp_unit = 'Mmol O' !'mmol O_2/m^3'
221	CASE ('mmol P/m^3') ; tmp_unit = 'Mmol P' !'mmol P/m^3'
222	CASE ('mmol Si/m^3') ; tmp_unit = 'Mmol S' !'mmol Si/m^3'
223	CASE ('umol/kg') ; tmp_unit = 'Mmol' !'umol/kg' = mmol/m^3
224	! CASE ('1/m') ; cycle
225	CASE DEFAULT
226	tmp_unit = TRIM(TRIM(tmp_unit)//TRIM('x 1e9 m^3'))
227	END SELECT
228	ENDIF
229
230	WRITE(numout,*) TRIM(TRIM('<field id="')//TRIM(tmp_name)//TRIM('" long_name="')// &
231	& TRIM(BGC_lev_name(jl))//TRIM('_')//TRIM(tmp_long_name)//TRIM(BGC_stat_name(js))// &
232	& TRIM('" unit="'//TRIM(tmp_unit) //'" />'))
233
234	END DO
235	END DO
236	END DO
237
238	WRITE(numout,*) ''
239	WRITE(numout,*) 'diaregmean BGC iodef.xml entries'
240	WRITE(numout,*) ''
241	DO js=1,6
242
243	DO jn=1,jp_fabm ! State loop
244
245	DO bgci=1,nBGC_output!
246	if (trim(ctrcnm(jn)) == TRIM(BGC_output_var(bgci))) CYCLE
247	ENDDO
248	DO jl=1,BGC_nlevs
249	! only print out area averages for ss, nb, diff, and depth averaged, and total values for volume integrated
250	IF ((jl .EQ. 5) .AND. (js .NE. 2)) CYCLE ! cycle if vol, and not tot.
251	IF ((jl .NE. 5) .AND. (js .NE. 1)) CYCLE ! cycle if other levels, and not ave.
252
253	tmp_name=TRIM(TRIM("reg_")//TRIM(BGC_lev_name(jl))//TRIM("_")//TRIM(ctrcnm(jn))// TRIM(BGC_stat_name(js)))
254	tmp_long_name = TRIM(ctrcln(jn))
255
256	WRITE(numout,*) TRIM(TRIM('<field field_ref="')//TRIM(tmp_name)//TRIM('"/>'))
257
258	END DO !level
259	END DO ! State loop
260	END DO !statistic
261	WRITE(numout,*) ''
262	DEALLOCATE( BGC_stat_name,BGC_lev_name)
263
264	ENDIF ! Control print
265
266	ENDIF !ln_diaregmean_bgc
267
268	#endif
269
270
271	IF (ln_diaregmean) THEN
272
273	! Open region mask for region means, and retrieve the size of the mask (number of levels)
274	CALL iom_open ( 'region_mask.nc', inum )
275	idmaskvar = iom_varid( inum, 'mask', kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE.)
276	nmasks = zdimsz(3)
277
278	! read in the region mask (which contains floating point numbers) into a temporary array of reals.
279	ALLOCATE( tmp_region_mask_real(jpi,jpj,nmasks), STAT= ierr )
280	IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate tmp_region_mask_real array' )
281
282	! Use jpdom_unknown to read in a n-layer mask.
283	tmp_region_mask_real(:,:,:) = 0
284	CALL iom_get( inum, jpdom_unknown, 'mask', tmp_region_mask_real(1:nlci,1:nlcj,1:nmasks), &
285	& kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nmasks /) )
286
287	CALL iom_close( inum )
288
289	!Convert the region mask of reals into one of integers.
290
291	ALLOCATE( region_mask(jpi,jpj,nmasks), STAT= ierr )
292	IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate region_mask array' )
293	region_mask(:,:,:) = 0
294	region_mask = int(tmp_region_mask_real(:,:,:))
295	DEALLOCATE( tmp_region_mask_real)
296
297
298	ALLOCATE( nreg_mat(nmasks), STAT= ierr )
299	IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate nreg_mat array' )
300
301	! work out the number of regions in each mask, asssuming land is 0, and the regions are consectively numbered,
302	! without missing any number, so the number of regions is the maximum number + 1 (for land). mpp_max across the
303	! processors to get the global maxima
304	DO maskno = 1,nmasks
305	tmpint = maxval(region_mask(:,:,maskno))
306	CALL mpp_max( tmpint )
307	nreg_mat(maskno) = tmpint + 1
308	END DO
309
310	IF(lwp) THEN
311	! if writing out as binary and text, open the files.
312	IF ( ln_diaregmean_bin ) THEN
313	! Open binary for region means
314	CALL ctl_opn( numdct_reg_bin ,'region_mean_timeseries.dat' , 'NEW', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
315	ENDIF
316
317	IF ( ln_diaregmean_ascii ) THEN
318	! Open text files for region means
319	CALL ctl_opn( numdct_reg_txt ,'region_mean_timeseries.txt' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
320	ENDIF
321	ENDIF
322	ENDIF
323
324	END SUBROUTINE dia_regmean_init
325
326	SUBROUTINE dia_calctmb_region_mean( pinfield,pouttmb )
327	!!---------------------------------------------------------------------
328	!! * ROUTINE dia_calctmb_region_mean *
329	!!
330	!! ** Purpose : Find the Top, Bottom and Top minus Bottom fields of water Column
331	!! : and depth average, and volume and mass intergated values.
332
333	!!
334	!! ** Method :
335	!! use mbathy to find surface, mid and bottom of model levels
336	!!
337	!! History :
338	!! 3.6 ! 08-14 (E. O'Dea) Routine based on dia_wri_foam
339	!!----------------------------------------------------------------------
340	!! * Modules used
341
342	! Routine to map 3d field to top, middle, bottom
343	IMPLICIT NONE
344
345
346	! Routine arguments
347	REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN ) :: pinfield ! Input 3d field and mask
348	REAL(wp), DIMENSION(jpi, jpj, 6 ), INTENT( OUT) :: pouttmb ! Output top, bottom and surface minus bed, zav, vol int, mass int
349
350	! Local variables
351	INTEGER :: ji,jj,jk ! Dummy loop indices
352
353	! Local Real
354	REAL(wp) :: zmdi ! set masked values
355	! for depth int
356	REAL(wp) :: tmpnumer,tmpnumer_mass,tmpdenom ,z_av_val,vol_int_val
357
358	zmdi=1.e+20 !missing data indicator for masking
359
360	!zmdi=0 !missing data indicator for masking
361
362	! Calculate top
363	pouttmb(:,:,1) = pinfield(:,:,1)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
364
365	! Calculate middle
366	!DO jj = 1,jpj
367	! DO ji = 1,jpi
368	! jk = max(1,mbathy(ji,jj)/2)
369	! pouttmb(ji,jj,2) = pinfield(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
370	! END DO
371	!END DO
372
373	! Calculate bottom, and top minus bottom
374	DO jj = 1,jpj
375	DO ji = 1,jpi
376	IF ( tmask(ji,jj,1) .EQ. 1) THEN ! if land
377
378	jk = max(1,mbathy(ji,jj) - 1)
379	pouttmb(ji,jj,2) = pinfield(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
380
381	pouttmb(ji,jj,3) = (pouttmb(ji,jj,1) - pouttmb(ji,jj,2))tmask(ji,jj,1) + zmdi(1.0-tmask(ji,jj,1))
382
383	!Depth and volume integral:
384	!---------------------------
385	!Vol int = Concentration * vol of grid box, summed over depth.
386	!Mass int = Concentration * vol of grid box * density of water, summed over depth.
387	!Depth Average = Vol int divided by * (vol of grid box summed over depth).
388
389	tmpnumer = 0.
390	tmpnumer_mass = 0.
391	tmpdenom = 0.
392	DO jk = 1,jpk
393	tmpnumer = tmpnumer + pinfield(ji,jj,jk)tmask(ji,jj,jk)e1t(ji,jj)e2t(ji,jj)e3t_n(ji,jj,jk)
394	tmpnumer_mass = tmpnumer_mass + pinfield(ji,jj,jk)tmask(ji,jj,jk)e1t(ji,jj)e2t(ji,jj)e3t_n(ji,jj,jk)*rhop(ji,jj,jk)
395	tmpdenom = tmpdenom + tmask(ji,jj,jk)e1t(ji,jj)e2t(ji,jj)*e3t_n(ji,jj,jk)
396	END DO
397	!z_av_val = tmpnumer/tmpdenom
398	!vol_int_val = tmpnumer
399	!mass_int_val = tmpnumer*density
400
401	pouttmb(ji,jj,4) = tmpnumer/tmpdenom ! depth averaged
402	pouttmb(ji,jj,5) = tmpnumer ! Vol integrated
403	pouttmb(ji,jj,6) = tmpnumer_mass ! Mass integrated (for heat and salt calcs)
404	ELSE
405	pouttmb(ji,jj,1) = zmdi
406	pouttmb(ji,jj,2) = zmdi
407	pouttmb(ji,jj,3) = zmdi
408	pouttmb(ji,jj,4) = zmdi
409	pouttmb(ji,jj,5) = zmdi
410	pouttmb(ji,jj,6) = zmdi
411	ENDIF
412	END DO
413	END DO
414
415	END SUBROUTINE dia_calctmb_region_mean
416
417
418	SUBROUTINE dia_regmean( kt )
419	!!----------------------------------------------------------------------
420	!! * ROUTINE dia_regmean *
421	!! ** Purpose : Produce regional mean diagnostics
422	!!
423	!! ** Method : calls dia_wri_region_mean to calculate and write the regional means for a number of variables,
424	!! (calling dia_calctmb_region_mean where necessary).
425	!!
426	!! Closes all text and binary files on last time step
427	!!
428	!!
429	!!
430	!!
431	!! History :
432	!! 3.6 ! 11-16 (J. Tinker)
433	!!
434	!!--------------------------------------------------------------------
435	REAL(wp), POINTER, DIMENSION(:,:,:) :: tmp1mat ! temporary array of 1's
436	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbT ! temporary T workspace
437	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbS ! temporary S workspace
438	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb1 ! temporary density workspace
439	REAL(wp) :: zmdi ! set masked values
440	INTEGER, INTENT( in ) :: kt ! ocean time-step index
441
442	REAL(wp) :: zdt ! temporary reals
443	INTEGER :: i_steps, ierr ! no of timesteps per hour, allocation error index
444	INTEGER :: maskno,jj,ji,jk,jm,nreg ! indices of mask, i and j, and number of regions
445
446	#if defined key_fabm
447	INTEGER :: jn ,tmp_dummy ! set masked values
448	REAL(wp) :: tmp_val ! tmp value, to allow min and max value clamping (not implemented)
449	INTEGER :: jl
450	CHARACTER (len=60) :: tmp_name_bgc_top,tmp_name_bgc_bot,tmp_name_bgc_dif, tmp_name_bgc_zav, tmp_name_bgc_vol
451	CHARACTER (len=60) :: tmp_output_filename
452	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbBGC ! temporary BGC workspace
453
454	LOGICAL :: verbose
455	verbose = .FALSE.
456	tmp_val = 0
457	#endif
458	zmdi=1.e+20 !missing data indicator for maskin
459
460	IF (ln_diaregmean) THEN
461	! If regional mean calculations required by namelist
462	! -----------------
463	! identify hourly time steps (not used)
464	zdt = rdt
465	IF( nacc == 1 ) zdt = rdtmin
466
467	IF( MOD( 3600,INT(zdt) ) == 0 ) THEN
468	i_steps = 3600/INT(zdt)
469	ELSE
470	CALL ctl_stop('STOP', 'dia_regmean: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
471	ENDIF
472
473	! Every time step, add physical, SBC, PEA, MLD terms to create hourly sums.
474	! Every hour, then hourly sums are divided by the number of timesteps in the hour to make hourly means
475	! These hourly mean values are then used to caluclate the regional means, and output with IOM.
476	#if defined key_fabm
477	! BGC values are not averaged up over the hour, but are output as hourly instantaneous values.
478	#endif
479
480
481	!Extract 2d fields from 3d T and S with dia_calctmb_region_mean
482	CALL wrk_alloc( jpi , jpj, 6 , zwtmbT )
483	CALL wrk_alloc( jpi , jpj, 6 , zwtmbS )
484	CALL wrk_alloc( jpi , jpj, 6 , zwtmb1 )
485
486	CALL dia_calctmb_region_mean( tsn(:,:,:,jp_tem),zwtmbT)
487	CALL dia_calctmb_region_mean( tsn(:,:,:,jp_sal),zwtmbS)
488
489	! To calc regional mean time series of int vol and mass, run region mean code on array of 1's...
490	! - then when multplying by volume, gives volume,
491	! - then when multplying by volume*density, gives mass
492
493	CALL wrk_alloc( jpi , jpj, jpk , tmp1mat )
494	DO jj = 1,jpj
495	DO ji = 1,jpi
496	DO jk = 1,jpk
497	tmp1mat(ji,jj,jk) = 1
498	END DO
499	END DO
500	END DO
501
502	CALL dia_calctmb_region_mean( tmp1mat,zwtmb1)
503	CALL wrk_dealloc( jpi , jpj, jpk , tmp1mat )
504
505	! Add 2d fields every time step to the hourly total.
506
507	tmp_field_mat(:,:,1) = tmp_field_mat(:,:,1) + (zwtmbT(:,:,1)*tmask(:,:,1)) !sst
508	tmp_field_mat(:,:,2) = tmp_field_mat(:,:,2) + (zwtmbT(:,:,2)*tmask(:,:,1)) !nbt
509	tmp_field_mat(:,:,3) = tmp_field_mat(:,:,3) + (zwtmbT(:,:,3)*tmask(:,:,1)) !dft
510
511	tmp_field_mat(:,:,4) = tmp_field_mat(:,:,4) + (zwtmbT(:,:,4)*tmask(:,:,1)) !zat
512	tmp_field_mat(:,:,5) = tmp_field_mat(:,:,5) + (zwtmbT(:,:,5)*tmask(:,:,1)) !vat
513	tmp_field_mat(:,:,6) = tmp_field_mat(:,:,6) + ((zwtmbT(:,:,6)tmask(:,:,1)4.2e3))! heat
514
515	tmp_field_mat(:,:,7) = tmp_field_mat(:,:,7) + (zwtmbS(:,:,1)*tmask(:,:,1)) !sss
516	tmp_field_mat(:,:,8) = tmp_field_mat(:,:,8) + (zwtmbS(:,:,2)*tmask(:,:,1)) !nbs
517	tmp_field_mat(:,:,9) = tmp_field_mat(:,:,9) + (zwtmbS(:,:,3)*tmask(:,:,1)) !dfs
518
519	tmp_field_mat(:,:,10) = tmp_field_mat(:,:,10) + (zwtmbS(:,:,4)*tmask(:,:,1)) !zas
520	tmp_field_mat(:,:,11) = tmp_field_mat(:,:,11) + (zwtmbS(:,:,5)*tmask(:,:,1)) !vas
521	tmp_field_mat(:,:,12) = tmp_field_mat(:,:,12) + (zwtmbS(:,:,6)*tmask(:,:,1)) !salt
522
523	tmp_field_mat(:,:,13) = tmp_field_mat(:,:,13) + (zwtmb1(:,:,5)*tmask(:,:,1))!vol
524	tmp_field_mat(:,:,14) = tmp_field_mat(:,:,14) + (zwtmb1(:,:,6)*tmask(:,:,1))!mass
525
526	tmp_field_mat(:,:,15) = tmp_field_mat(:,:,15) + (sshn(:,:)*tmask(:,:,1)) !ssh
527
528	CALL wrk_dealloc( jpi , jpj, 6 , zwtmbT )
529	CALL wrk_dealloc( jpi , jpj, 6 , zwtmbS )
530	CALL wrk_dealloc( jpi , jpj, 6 , zwtmb1 )
531
532	IF( ln_diaregmean_karamld ) THEN
533	tmp_field_mat(:,:,16) = tmp_field_mat(:,:,16) + (hmld_zint(:,:)*tmask(:,:,1)) !mldkara
534	ENDIF
535
536	IF( ln_diaregmean_pea ) THEN
537	tmp_field_mat(:,:,17) = tmp_field_mat(:,:,17) + (pea(:,:)*tmask(:,:,1)) !pea
538	tmp_field_mat(:,:,18) = tmp_field_mat(:,:,18) + (peat(:,:)*tmask(:,:,1)) !peat
539	tmp_field_mat(:,:,19) = tmp_field_mat(:,:,19) + (peas(:,:)*tmask(:,:,1)) !peas
540	ENDIF
541
542	IF( ln_diaregmean_diaar5 ) THEN
543	tmp_field_AR5_mat(:,:,1) = tmp_field_AR5_mat(:,:,1) + (sshsteric_mat(:,:)*tmask(:,:,1))
544	tmp_field_AR5_mat(:,:,2) = tmp_field_AR5_mat(:,:,2) + (sshthster_mat(:,:)*tmask(:,:,1))
545	tmp_field_AR5_mat(:,:,3) = tmp_field_AR5_mat(:,:,3) + (sshhlster_mat(:,:)*tmask(:,:,1))
546	tmp_field_AR5_mat(:,:,4) = tmp_field_AR5_mat(:,:,4) + (zbotpres_mat(:,:)*tmask(:,:,1))
547	ENDIF
548
549	IF( ln_diaregmean_diasbc ) THEN
550	tmp_field_SBC_mat(:,:,1) = tmp_field_SBC_mat(:,:,1) + ((qsr + qns)*tmask(:,:,1))
551	tmp_field_SBC_mat(:,:,2) = tmp_field_SBC_mat(:,:,2) + (qsr*tmask(:,:,1))
552	tmp_field_SBC_mat(:,:,3) = tmp_field_SBC_mat(:,:,3) + (qns*tmask(:,:,1))
553	tmp_field_SBC_mat(:,:,4) = tmp_field_SBC_mat(:,:,4) + (emp*tmask(:,:,1))
554	tmp_field_SBC_mat(:,:,5) = tmp_field_SBC_mat(:,:,5) + (wndm*tmask(:,:,1))
555	tmp_field_SBC_mat(:,:,6) = tmp_field_SBC_mat(:,:,6) + (pressnow*tmask(:,:,1))
556	tmp_field_SBC_mat(:,:,7) = tmp_field_SBC_mat(:,:,7) + (rnf*tmask(:,:,1))
557	ENDIF
558
559
560
561	tmp_field_cnt = tmp_field_cnt + 1
562
563	! On the hour, calculate hourly means from the hourly total,and process the regional means.
564
565	IF ( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 ) THEN
566
567
568	CALL dia_wri_region_mean(kt, "sst" , tmp_field_mat(:,:,1)/real(tmp_field_cnt,wp))
569	CALL dia_wri_region_mean(kt, "nbt" , tmp_field_mat(:,:,2)/real(tmp_field_cnt,wp))
570	CALL dia_wri_region_mean(kt, "dft" , tmp_field_mat(:,:,3)/real(tmp_field_cnt,wp))
571
572	CALL dia_wri_region_mean(kt, "zat" , tmp_field_mat(:,:,4)/real(tmp_field_cnt,wp))
573	CALL dia_wri_region_mean(kt, "vat" , tmp_field_mat(:,:,5)/real(tmp_field_cnt,wp))
574	CALL dia_wri_region_mean(kt, "heat" , tmp_field_mat(:,:,6)/real(tmp_field_cnt,wp)/1e12)
575
576	CALL dia_wri_region_mean(kt, "sss" , tmp_field_mat(:,:,7)/real(tmp_field_cnt,wp))
577	CALL dia_wri_region_mean(kt, "nbs" , tmp_field_mat(:,:,8)/real(tmp_field_cnt,wp))
578	CALL dia_wri_region_mean(kt, "dfs" , tmp_field_mat(:,:,9)/real(tmp_field_cnt,wp))
579
580	CALL dia_wri_region_mean(kt, "zas" , tmp_field_mat(:,:,10)/real(tmp_field_cnt,wp))
581	CALL dia_wri_region_mean(kt, "vas" , tmp_field_mat(:,:,11)/real(tmp_field_cnt,wp))
582	CALL dia_wri_region_mean(kt, "salt" , tmp_field_mat(:,:,12)/real(tmp_field_cnt,wp)/1e12)
583
584	CALL dia_wri_region_mean(kt, "vol" , tmp_field_mat(:,:,13)/real(tmp_field_cnt,wp))
585	CALL dia_wri_region_mean(kt, "mass" , tmp_field_mat(:,:,14)/real(tmp_field_cnt,wp))
586
587	CALL dia_wri_region_mean(kt, "ssh" , tmp_field_mat(:,:,15)/real(tmp_field_cnt,wp))
588
589
590	IF( ln_diaregmean_karamld ) THEN
591	CALL dia_wri_region_mean(kt, "mldkara" , tmp_field_mat(:,:,16)/real(tmp_field_cnt,wp)) ! tm
592	ENDIF
593
594	IF( ln_diaregmean_pea ) THEN
595	CALL dia_wri_region_mean(kt, "pea" , tmp_field_mat(:,:,17)/real(tmp_field_cnt,wp))
596	CALL dia_wri_region_mean(kt, "peat" , tmp_field_mat(:,:,18)/real(tmp_field_cnt,wp))
597	CALL dia_wri_region_mean(kt, "peas" , tmp_field_mat(:,:,19)/real(tmp_field_cnt,wp)) ! tmb
598	ENDIF
599
600	tmp_field_mat(:,:,:) = 0.
601
602	IF( ln_diaregmean_diaar5 ) THEN
603
604	CALL dia_wri_region_mean(kt, "ssh_steric" , tmp_field_AR5_mat(:,:,1)/real(tmp_field_cnt,wp))
605	CALL dia_wri_region_mean(kt, "ssh_thermosteric", tmp_field_AR5_mat(:,:,2)/real(tmp_field_cnt,wp))
606	CALL dia_wri_region_mean(kt, "ssh_halosteric" , tmp_field_AR5_mat(:,:,3)/real(tmp_field_cnt,wp))
607	CALL dia_wri_region_mean(kt, "bot_pres" , tmp_field_AR5_mat(:,:,4)/real(tmp_field_cnt,wp))
608	tmp_field_AR5_mat(:,:,:) = 0.
609	ENDIF
610
611	IF( ln_diaregmean_diasbc ) THEN
612
613	CALL dia_wri_region_mean(kt, "qt" , tmp_field_SBC_mat(:,:,1)/real(tmp_field_cnt,wp))
614	CALL dia_wri_region_mean(kt, "qsr" , tmp_field_SBC_mat(:,:,2)/real(tmp_field_cnt,wp))
615	CALL dia_wri_region_mean(kt, "qns" , tmp_field_SBC_mat(:,:,3)/real(tmp_field_cnt,wp))
616	CALL dia_wri_region_mean(kt, "emp" , tmp_field_SBC_mat(:,:,4)/real(tmp_field_cnt,wp))
617	CALL dia_wri_region_mean(kt, "wspd" , tmp_field_SBC_mat(:,:,5)/real(tmp_field_cnt,wp))
618	CALL dia_wri_region_mean(kt, "mslp" , tmp_field_SBC_mat(:,:,6)/real(tmp_field_cnt,wp))
619	CALL dia_wri_region_mean(kt, "rnf" , tmp_field_SBC_mat(:,:,7)/real(tmp_field_cnt,wp))
620	tmp_field_SBC_mat(:,:,:) = 0.
621	ENDIF
622
623	#if defined key_fabm
624	!ADD Biogeochemistry
625
626	IF( ln_diaregmean_bgc ) THEN !ln_diaregmean_bgc
627
628	! Loop through 3d BGC tracers
629	DO jn=1,jp_fabm ! State loop
630
631	! get variable name for different levels
632	tmp_name_bgc_top=TRIM(TRIM("top_")//TRIM(ctrcnm(jn)))
633	tmp_name_bgc_bot=TRIM(TRIM("bot_")//TRIM(ctrcnm(jn)))
634	tmp_name_bgc_dif=TRIM(TRIM("dif_")//TRIM(ctrcnm(jn)))
635	tmp_name_bgc_zav=TRIM(TRIM("zav_")//TRIM(ctrcnm(jn)))
636	tmp_name_bgc_vol=TRIM(TRIM("vol_")//TRIM(ctrcnm(jn)))
637
638	! print out names if verbose
639	IF(verbose .AND. lwp) THEN
640	WRITE(numout,*)
641	WRITE(numout,*) 'dia_regmean tmp_name_bgc_top : ',TRIM(tmp_name_bgc_top)
642	WRITE(numout,*) 'dia_regmean tmp_name_bgc_bot : ',TRIM(tmp_name_bgc_bot)
643	WRITE(numout,*) 'dia_regmean tmp_name_bgc_dif : ',TRIM(tmp_name_bgc_dif)
644	WRITE(numout,*) 'dia_regmean tmp_name_bgc_zav : ',TRIM(tmp_name_bgc_zav)
645	WRITE(numout,*) 'dia_regmean tmp_name_bgc_vol : ',TRIM(tmp_name_bgc_vol)
646	CALL FLUSH(numout)
647
648	ENDIF
649
650	!Allocate working array, and get surface, bed etc fields.
651	CALL wrk_alloc( jpi , jpj, 6 , zwtmbBGC )
652	CALL dia_calctmb_region_mean( trn(:,:,:,jn),zwtmbBGC )
653
654
655	!Print out 2d fields to ascii text files to check values if verbose. (24MB per time step, per BGC variable)
656	IF (verbose) THEN
657
658	WRITE (tmp_output_filename, "(A4,I3.3,A1,I6.6,A1,I3.3,A4)") "bgc_",jn,"_",kt,"_",narea,".txt"
659	WRITE (,) tmp_output_filename
660	OPEN(UNIT=74,FILE=TRIM(tmp_output_filename))
661
662	DO ji = 1,jpi
663	DO jj = 1,jpj
664	WRITE(74,FMT="(I4,I4,F3,F25.5,F25.5,F25.5,F25.5,F25.5)") nimpp+ji, njmpp+jj,tmask(ji,jj,1),&
665	& zwtmbBGC(ji,jj,1),zwtmbBGC(ji,jj,2),zwtmbBGC(ji,jj,3),zwtmbBGC(ji,jj,4),zwtmbBGC(ji,jj,5)/1e9
666	END DO
667	END DO
668	CLOSE(74)
669	ENDIF
670
671	! Do region means
672	CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_top) , zwtmbBGC(:,:,1))
673	CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_bot) , zwtmbBGC(:,:,2))
674	CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_dif) , zwtmbBGC(:,:,3))
675	CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_zav) , zwtmbBGC(:,:,4))
676	CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_vol) , zwtmbBGC(:,:,5)/1e9)
677
678
679	!Deallocate working array
680	CALL wrk_dealloc( jpi , jpj, 6 , zwtmbBGC )
681	ENDDO ! State loop
682	ENDIF !ln_diaregmean_bgc
683
684	#endif
685
686	tmp_field_cnt = 0
687
688	ENDIF ! ( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 )
689
690
691	! If on the last time step, close binary and ascii files.
692	IF( kt == nitend ) THEN
693	IF(lwp) THEN
694	IF ( ln_diaregmean_bin ) THEN
695	!Closing binary files for regional mean time series.
696	CLOSE(numdct_reg_bin)
697	ENDIF
698	IF ( ln_diaregmean_ascii ) THEN
699	!Closing text files for regional mean time series.
700	CLOSE(numdct_reg_txt)
701	ENDIF
702
703	DEALLOCATE( region_mask, nreg_mat, tmp_field_mat)
704	IF( ln_diaregmean_diaar5 ) DEALLOCATE( tmp_field_AR5_mat)
705	IF( ln_diaregmean_diasbc ) DEALLOCATE( tmp_field_SBC_mat)
706	ENDIF
707	ENDIF
708
709
710	ELSE
711	CALL ctl_warn('dia_regmean: regmean diagnostic is set to false you should not have seen this')
712	ENDIF
713
714	END SUBROUTINE dia_regmean
715
716
717	SUBROUTINE dia_wri_region_mean(kt, tmp_name, infield )
718	!!---------------------------------------------------------------------
719	!! * ROUTINE dia_tmb *
720	!!
721	!! ** Purpose : Calculate and write region mean time series for 2d arrays
722	!!
723	!! ** Method :
724	!! use
725	!!
726	!! History :
727	!! ?? ! 15/10/2015 (JTinker) Routine taken from old dia_wri_foam
728	!!----------------------------------------------------------------------
729	!! * Modules used
730	!use lib_mpp
731	!use lib_fortr
732	IMPLICIT NONE
733
734	INTEGER, INTENT(in) :: kt
735	CHARACTER (len=*) , INTENT(IN ) :: tmp_name
736	REAL(wp), DIMENSION(jpi, jpj), INTENT(IN ) :: infield ! Input 3d field and mask
737
738	! Local variables
739	INTEGER, DIMENSION(jpi, jpj) :: internal_region_mask ! Input 3d field and mask
740	REAL(wp), DIMENSION(jpi, jpj) :: internal_infield ! Internal data field
741	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zrmet_ave,zrmet_tot,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id ,zrmet_min,zrmet_max
742	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zrmet_out
743	REAL(wp), ALLOCATABLE, DIMENSION(:) :: ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,reg_id_mat,mask_id_mat ,min_mat,max_mat !: region_mask
744
745	REAL(wp) :: zmdi, zrmet_val ! set masked values
746	INTEGER :: maskno,nreg ! ocean time-step indexocean time step
747	INTEGER :: ji,jj,jk,ind,jm ! Dummy loop indices
748	INTEGER :: reg_ind_cnt ! Dummy loop indices
749
750	INTEGER :: ierr
751	REAL(wp) :: tmpreal
752	CHARACTER(LEN=180) :: FormatString,nreg_string,tmp_name_iom
753	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: dummy_zrmet
754	LOGICAL :: verbose
755	verbose = .false.
756
757
758	zmdi=1.e+20 !missing data indicator for maskin
759
760	!Allocate output arrays for iomput, set to zmdi, and set a region counter = 1
761	ALLOCATE( zrmet_ave(n_regions_output), STAT= ierr )
762	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_ave array' )
763	ALLOCATE( zrmet_tot(n_regions_output), STAT= ierr )
764	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_tot array' )
765	ALLOCATE( zrmet_var(n_regions_output), STAT= ierr )
766	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_var array' )
767	ALLOCATE( zrmet_cnt(n_regions_output), STAT= ierr )
768	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_cnt array' )
769	ALLOCATE( zrmet_mask_id(n_regions_output), STAT= ierr )
770	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_mask_id array' )
771	ALLOCATE( zrmet_reg_id(n_regions_output), STAT= ierr )
772	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
773
774
775	ALLOCATE( zrmet_min(n_regions_output), STAT= ierr )
776	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_min array' )
777	ALLOCATE( zrmet_max(n_regions_output), STAT= ierr )
778	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_max array' )
779
780	ALLOCATE( zrmet_out(jpi,jpj,n_regions_output), STAT= ierr )
781	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
782
783
784
785	IF(lwp .AND. verbose) THEN
786	WRITE(numout,*)
787	WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//';'
788	WRITE(numout,*)
789	ENDIF
790
791	DO ji = 1,jpi
792	DO jj = 1,jpj
793	internal_infield(ji,jj) = infield(ji,jj)
794	END DO
795	END DO
796
797	! Check for NANS # JT 03/09/2018
798	DO ji = 1,jpi
799	DO jj = 1,jpj
800	IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
801	IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
802	WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//' Nan at (kt,i,j): ',kt,ji - (-jpizoom+1-nimpp+1),jj - (-jpjzoom+1-njmpp+1)
803	internal_infield(ji,jj) = 0.
804	ENDIF
805	ELSE
806	IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
807	WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//' Masked Nan at (kt,i,j): ',kt,ji - (-jpizoom+1-nimpp+1),jj - (-jpjzoom+1-njmpp+1)
808	internal_infield(ji,jj) = 0.
809	ENDIF
810	ENDIF
811	END DO
812	END DO
813
814
815	zrmet_ave(:) = zmdi
816	zrmet_tot(:) = zmdi
817	zrmet_var(:) = zmdi
818	zrmet_cnt(:) = zmdi
819	zrmet_mask_id(:) = zmdi
820	zrmet_reg_id(:) = zmdi
821
822	zrmet_min(:) = zmdi
823	zrmet_max(:) = zmdi
824	reg_ind_cnt = 1
825
826
827	! loop though the masks
828	DO maskno = 1,nmasks
829	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin mask loops: ',maskno
830
831
832	! For each mask, get the number of regions (nreg), and a local copy of the region.
833	nreg = nreg_mat(maskno)
834	internal_region_mask = region_mask(:,:,maskno)
835
836	! allocate temporary stat arrays, and set to zero
837	ALLOCATE( ave_mat(nreg), STAT= ierr )
838	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ave_mat array' )
839	ALLOCATE( tot_mat(nreg), STAT= ierr )
840	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate tot_mat array' )
841	ALLOCATE( num_mat(nreg), STAT= ierr )
842	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate num_mat array' )
843	ALLOCATE( var_mat(nreg), STAT= ierr )
844	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate var_mat array' )
845	ALLOCATE( ssq_mat(nreg), STAT= ierr )
846	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ssq_mat array' )
847	ALLOCATE( cnt_mat(nreg), STAT= ierr )
848	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate cnt_mat array' )
849
850	ALLOCATE( min_mat(nreg), STAT= ierr )
851	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate min_mat array' )
852	ALLOCATE( max_mat(nreg), STAT= ierr )
853	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate max_mat array' )
854
855	ALLOCATE( reg_id_mat(nreg), STAT= ierr )
856	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate reg_id_mat array' )
857	ALLOCATE( mask_id_mat(nreg), STAT= ierr )
858	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate mask_id_mat array' )
859
860
861
862	ave_mat(:) = 0.
863	tot_mat(:) = 0.
864	num_mat(:) = 0.
865	var_mat(:) = 0.
866	cnt_mat(:) = 0.
867	ssq_mat(:) = 0.
868
869	min_mat(:) = zmdi
870	max_mat(:) = -zmdi
871	reg_id_mat(:) = 0.
872	mask_id_mat(:) = 0.
873
874	! loop though the array. for each sea grid box where tmask == 1),
875	! read which region the grid box is in, add the value of the gridbox (and its square)
876	! to the total for that region, and then increment the counter for that region.
877	!CALL cpu_time(start_reg_mean_loop)
878	!WRITE(numout,*) kt,start_reg_mean_loop
879	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin spatial loops: '
880	DO ji = nldi,nlei
881	DO jj = nldj,nlej
882	IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
883	ind = internal_region_mask(ji,jj)+1
884	tot_mat(ind) = tot_mat(ind) + (internal_infield(ji,jj))
885	ssq_mat(ind) = ssq_mat(ind) + ( internal_infield(ji,jj) * internal_infield(ji,jj))
886	cnt_mat(ind) = cnt_mat(ind) + 1.
887
888	min_mat(ind) = min(min_mat(ind),internal_infield(ji,jj))
889	max_mat(ind) = max(max_mat(ind),internal_infield(ji,jj))
890	ENDIF
891	END DO
892	END DO
893	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finish spatial loops: '
894	! sum the totals, the counts, and the squares across the processors
895	CALL mpp_sum( tot_mat,nreg )
896	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum 1'
897	CALL mpp_sum( ssq_mat,nreg )
898	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum 2'
899	CALL mpp_sum( cnt_mat,nreg )
900	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum 2'
901
902	CALL mpp_min( min_mat,nreg )
903	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_min'
904	CALL mpp_max( max_mat,nreg )
905	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_max'
906
907
908	!calculate the mean and variance from the total, sum of squares and the count.
909
910	ave_mat = tot_mat(:)/cnt_mat(:)
911	var_mat = ssq_mat(:)/cnt_mat(:) - (ave_mat(:)*ave_mat(:))
912
913
914	!mask array of mask and region number.
915	DO jj = 1,nreg
916	reg_id_mat(jj) = real(jj-1)
917	mask_id_mat(jj) = real(maskno)
918	END DO
919
920
921	!write text and binary, and note region statistics for current mask for later iom_put
922	IF( lwp ) THEN
923
924	!Write out ascii and binary if requred
925	IF ( ln_diaregmean_bin ) THEN
926	!Writing out regional mean time series to binary files
927	WRITE(numdct_reg_bin) tmp_name,kt,maskno,n_regions_output
928	WRITE(numdct_reg_bin) ave_mat
929	WRITE(numdct_reg_bin) tot_mat
930	WRITE(numdct_reg_bin) var_mat
931	WRITE(numdct_reg_bin) ssq_mat
932	WRITE(numdct_reg_bin) cnt_mat
933	WRITE(numdct_reg_bin) min_mat
934	WRITE(numdct_reg_bin) max_mat
935	ENDIF
936
937	IF ( ln_diaregmean_ascii ) THEN
938	!Writing out regional mean time series to text files
939
940	WRITE(nreg_string, "(I5)") nreg
941	FormatString = "(A30,"//trim(nreg_string)//"F25.3)"
942	WRITE(numdct_reg_txt, FMT="(A30,I6,I6)") tmp_name,kt,maskno
943	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ave_mat:", ave_mat
944	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"tot_mat:", tot_mat
945	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"var_mat:", var_mat
946	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ssq_mat:", ssq_mat
947	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"cnt_mat:", cnt_mat
948	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"min_mat:", min_mat
949	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"max_mat:", max_mat
950	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"reg_mat:", reg_id_mat
951	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"msk_mat:", mask_id_mat
952
953	ENDIF
954
955	DO jm = 1,nreg
956	zrmet_ave( reg_ind_cnt) = ave_mat(jm)
957	zrmet_tot( reg_ind_cnt) = tot_mat(jm)
958	zrmet_var( reg_ind_cnt) = var_mat(jm)
959	zrmet_cnt( reg_ind_cnt) = cnt_mat(jm)
960	zrmet_min( reg_ind_cnt) = min_mat(jm)
961	zrmet_max( reg_ind_cnt) = max_mat(jm)
962	zrmet_reg_id( reg_ind_cnt) = reg_id_mat(jm)
963	zrmet_mask_id(reg_ind_cnt) = mask_id_mat(jm)
964
965	reg_ind_cnt = reg_ind_cnt + 1
966	END DO
967
968	ENDIF
969
970	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean about to deallocated arrays for ',kt,maskno
971	DEALLOCATE(ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,min_mat,max_mat,reg_id_mat,mask_id_mat)
972
973	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean deallocated arrays for ',kt,maskno
974	IF(lwp .AND. ln_diaregmean_ascii) CALL FLUSH(numdct_reg_txt)
975	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean flushed region mean text for ',kt,maskno
976	END DO
977
978	IF(lwp .AND. verbose) THEN ! Control print
979	WRITE(numout,*) 'dia_regmean ready to start iom_put'
980	CALL FLUSH(numout)
981	ENDIF
982
983	!With current field_def.xml and iodef.xml, these fields must be output, so set to dummy values if not required.
984
985	IF ( ln_diaregmean_nc ) THEN
986
987	zrmet_out(:,:,:) = 0
988	zrmet_val = 0
989	tmp_name_iom = ''
990
991	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean ready to start iom_put: ',trim(tmp_name)
992
993
994	DO jm = 1,n_regions_output
995	zrmet_val = zrmet_ave(jm)
996	! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
997	! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
998	if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
999	zrmet_out(:,:,jm) = zrmet_val
1000	END DO
1001	tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_ave'))
1002	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
1003	CALL iom_put(trim(tmp_name_iom), zrmet_out )
1004	zrmet_out(:,:,:) = 0
1005	zrmet_val = 0
1006	tmp_name_iom = ''
1007
1008	DO jm = 1,n_regions_output
1009	zrmet_val = zrmet_tot(jm)
1010	! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1011	! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1012	if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1013	zrmet_out(:,:,jm) = zrmet_val
1014	END DO
1015	tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_tot'))
1016	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
1017	CALL iom_put( trim(tmp_name_iom), zrmet_out )
1018	zrmet_out(:,:,:) = 0
1019	zrmet_val = 0
1020	tmp_name_iom = ''
1021
1022	DO jm = 1,n_regions_output
1023	zrmet_val = zrmet_var(jm)
1024	! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1025	! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1026	if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1027	zrmet_out(:,:,jm) = zrmet_val
1028	END DO
1029	tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_var'))
1030	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
1031	CALL iom_put( trim(tmp_name_iom), zrmet_out )
1032	zrmet_out(:,:,:) = 0
1033	zrmet_val = 0
1034	tmp_name_iom = ''
1035
1036	DO jm = 1,n_regions_output
1037	zrmet_val = zrmet_cnt(jm)
1038	! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1039	! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1040	if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1041	zrmet_out(:,:,jm) = zrmet_val
1042	END DO
1043	tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))
1044	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
1045	CALL iom_put( trim(tmp_name_iom), zrmet_out )
1046	zrmet_out(:,:,:) = 0
1047	zrmet_val = 0
1048	tmp_name_iom = ''
1049
1050	DO jm = 1,n_regions_output
1051	zrmet_val = zrmet_reg_id(jm)
1052	! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1053	! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1054	if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1055	zrmet_out(:,:,jm) = zrmet_val
1056	END DO
1057	tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_reg_id'))
1058	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
1059	CALL iom_put( trim(tmp_name_iom), zrmet_out )
1060	zrmet_out(:,:,:) = 0
1061	zrmet_val = 0
1062	tmp_name_iom = ''
1063
1064	DO jm = 1,n_regions_output
1065	zrmet_val = zrmet_mask_id(jm)
1066	! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1067	! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1068	if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1069	zrmet_out(:,:,jm) = zrmet_val
1070	END DO
1071	tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_mask_id'))
1072	IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
1073	CALL iom_put( trim(tmp_name_iom), zrmet_out )
1074	zrmet_out(:,:,:) = 0
1075	zrmet_val = 0
1076	tmp_name_iom = ''
1077	ELSE
1078
1079	ALLOCATE( dummy_zrmet(jpi,jpj,n_regions_output), STAT= ierr )
1080	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate dummy_zrmet array' )
1081
1082	DO jm = 1,n_regions_output
1083	dummy_zrmet(:,:,jm) = real(jm,wp)
1084	END DO
1085
1086	DO jm = 1,9
1087	CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_ave')), dummy_zrmet )
1088	CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_tot')), dummy_zrmet )
1089	CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_var')), dummy_zrmet )
1090	CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_cnt')), dummy_zrmet )
1091	CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_reg_id')), dummy_zrmet )
1092	CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_mask_id')), dummy_zrmet )
1093	END DO
1094
1095	DEALLOCATE( dummy_zrmet)
1096	ENDIF
1097
1098	DEALLOCATE(zrmet_ave,zrmet_tot,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id,zrmet_min,zrmet_max,zrmet_out)
1099
1100	IF(lwp .AND. verbose) THEN ! Control print
1101	WRITE(numout,*)
1102	WRITE(numout,*) 'dia_wri_region_mean finished for ', trim(tmp_name)
1103	WRITE(numout,*)
1104	CALL FLUSH(numout)
1105	ENDIF
1106
1107	END SUBROUTINE dia_wri_region_mean
1108
1109
1110	!!======================================================================
1111	END MODULE diaregmean

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

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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis3/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90

Download in other formats: