Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diaregmean.F90 @ 10492

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

Download in other formats: