Context Navigation

diaregmean.F90 @ 10489

Last change on this file since 10489 was 10489, checked in by hadjt, 5 years ago
Updated DIA/diaregmean.F90 to fix bug caused by double counting grid boxes where processors overlap
File size: 31.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 diatmb ! Top,middle,bottom output
14	USE diapea ! Top,middle,bottom output
15	USE zdfmxl ! MLD
16	USE sbc_oce
17	#if defined key_diaar5
18	USE diaar5
19	#endif
20	IMPLICIT NONE
21	PRIVATE
22
23	LOGICAL , PUBLIC :: ln_diaregmean ! region mean calculation
24	PUBLIC dia_regmean_init ! routine called by nemogcm.F90
25	PUBLIC dia_regmean ! routine called by diawri.F90
26
27
28
29
30
31	LOGICAL :: ln_diaregmean_ascii ! region mean calculation ascii output
32	LOGICAL :: ln_diaregmean_bin ! region mean calculation binary output
33	LOGICAL :: ln_diaregmean_nc ! region mean calculation netcdf output
34	LOGICAL :: ln_diaregmean_diaar5 ! region mean calculation including AR5 SLR terms
35	LOGICAL :: ln_diaregmean_diasbc ! region mean calculation including Surface BC
36	LOGICAL :: ln_diaregmean_karamld ! region mean calculation including kara mld terms
37	LOGICAL :: ln_diaregmean_pea ! region mean calculation including pea terms
38
39	REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: tmp_region_mask_real ! tempory region_mask of reals
40	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: region_mask ! region_mask matrix
41	INTEGER :: nmasks ! Number of mask files in region_mask.nc file -
42	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: nreg_mat ! Number of regions in each mask
43
44	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_field_mat !: temporary region_mask
45	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_field_AR5_mat !: temporary region_mask
46	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_field_SBC_mat !: temporary region_mask
47	INTEGER :: tmp_field_cnt ! tmp_field_cnt integer
48	!!----------------------------------------------------------------------
49	!! NEMO/OPA 3.6 , NEMO Consortium (2014)
50	!! $Id$
51	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
52	!!----------------------------------------------------------------------
53	CONTAINS
54
55	SUBROUTINE dia_regmean_init
56	!!---------------------------------------------------------------------------
57	!! * ROUTINE dia_regmean_init *
58	!!
59	!! ** Purpose: Initialization of region mask namelist
60	!!
61	!! ** Method : Read namelist
62	!! History
63	!! 3.6 ! 11-16 (J Tinker) Routine to initialize dia_regmean
64	!!---------------------------------------------------------------------------
65	!!
66	INTEGER :: ios ! Local integer output status for namelist read
67	INTEGER :: inum ! temporary logical unit ! copied from DOM/domzgr.F90
68	INTEGER :: ierr ! error integer for IOM_get
69	INTEGER :: idmaskvar ! output of iom_varid
70	INTEGER :: maskno ! counter for number of masks
71	INTEGER :: jj,ji ! i and j index
72	INTEGER :: tmpint ! temporary integer
73	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tmpregion !: temporary region_mask
74	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
75	INTEGER :: zndims ! number of dimensions in an array (i.e. 3, )
76	!
77	NAMELIST/nam_diaregmean/ ln_diaregmean,ln_diaregmean_ascii,ln_diaregmean_bin,ln_diaregmean_nc,&
78	& ln_diaregmean_karamld, ln_diaregmean_pea,ln_diaregmean_diaar5,ln_diaregmean_diasbc
79
80
81	! read in Namelist.
82	!!----------------------------------------------------------------------
83	!
84	REWIND ( numnam_ref ) ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
85	READ ( numnam_ref, nam_diaregmean, IOSTAT=ios, ERR= 901 )
86	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in reference namelist', lwp )
87
88	REWIND( numnam_cfg ) ! Namelist nam_diatmb in configuration namelist TMB diagnostics
89	READ ( numnam_cfg, nam_diaregmean, IOSTAT = ios, ERR = 902 )
90	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in configuration namelist', lwp )
91	IF(lwm) WRITE ( numond, nam_diaregmean )
92
93	IF(lwp) THEN ! Control print
94	WRITE(numout,*)
95	WRITE(numout,*) 'dia_regmean_init : Output regional mean Diagnostics'
96	WRITE(numout,*) '~~~~~~~~~~~~'
97	WRITE(numout,*) 'Namelist nam_regmean : set regmeanoutputs '
98	WRITE(numout,*) 'Switch for regmean diagnostics (T) or not (F) ln_diaregmean = ', ln_diaregmean
99	WRITE(numout,*) 'Switch for regmean ascii output (T) or not (F) ln_diaregmean_ascii = ', ln_diaregmean_ascii
100	WRITE(numout,*) 'Switch for regmean binary output (T) or not (F) ln_diaregmean_bin = ', ln_diaregmean_bin
101	WRITE(numout,*) 'Switch for regmean netcdf output (T) or not (F) ln_diaregmean_nc = ', ln_diaregmean_nc
102	WRITE(numout,*) 'Switch for regmean kara mld terms (T) or not (F) ln_diaregmean_karamld = ', ln_diaregmean_karamld
103	WRITE(numout,*) 'Switch for regmean PEA terms (T) or not (F) ln_diaregmean_pea = ', ln_diaregmean_pea
104	WRITE(numout,*) 'Switch for regmean AR5 SLR terms (T) or not (F) ln_diaregmean_diaar5 = ', ln_diaregmean_diaar5
105	WRITE(numout,*) 'Switch for regmean Surface forcing terms (T) or not (F) ln_diaregmean_diasbc = ', ln_diaregmean_diasbc
106	ENDIF
107
108
109	!ALLOCATE( tmp_field_mat(jpi,jpj,7), STAT= ierr ) !SS/NB/DT T/S, SSH, MLD, PEA, PEAT, PEAS
110	ALLOCATE( tmp_field_mat(jpi,jpj,11), STAT= ierr ) !SS/NB/DT T/S, SSH, MLD, PEA, PEAT, PEAS
111	IF( ierr /= 0 ) CALL ctl_stop( 'tmp_field_mat: failed to allocate tmp_region_mask_real array' )
112	tmp_field_mat(:,:,:) = 0.
113	tmp_field_cnt = 0
114
115	IF(ln_diaregmean_diaar5) THEN
116	ALLOCATE( tmp_field_AR5_mat(jpi,jpj,4), STAT= ierr ) !SS/NB/DT T/S, SSH, MLD, PEA, PEAT, PEAS
117	IF( ierr /= 0 ) CALL ctl_stop( 'tmp_field_AR5_mat: failed to allocate tmp_region_mask_real array' )
118	tmp_field_AR5_mat(:,:,:) = 0.
119	ENDIF
120
121	IF(ln_diaregmean_diasbc) THEN
122	ALLOCATE( tmp_field_SBC_mat(jpi,jpj,7), STAT= ierr ) !SS/NB/DT T/S, SSH, MLD, PEA, PEAT, PEAS
123	IF( ierr /= 0 ) CALL ctl_stop( 'tmp_field_SBC_mat: failed to allocate tmp_region_mask_real array' )
124	tmp_field_SBC_mat(:,:,:) = 0.
125	ENDIF
126
127	IF (ln_diaregmean) THEN
128
129	! Open region mask for region means, and retrieve the size of the mask (number of levels)
130	CALL iom_open ( 'region_mask.nc', inum )
131	idmaskvar = iom_varid( inum, 'mask', kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE.)
132	nmasks = zdimsz(3)
133
134	! read in the region mask (which contains floating point numbers) into a temporary array of reals.
135	ALLOCATE( tmp_region_mask_real(jpi,jpj,nmasks), STAT= ierr )
136	IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate tmp_region_mask_real array' )
137
138	! Use jpdom_unknown to read in a n layer mask.
139	tmp_region_mask_real(:,:,:) = 0
140	CALL iom_get( inum, jpdom_unknown, 'mask', tmp_region_mask_real(1:nlci,1:nlcj,1:nmasks), &
141	& kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nmasks /) )
142
143	CALL iom_close( inum )
144
145	!Convert the region mask of reals into one of integers.
146
147	ALLOCATE( region_mask(jpi,jpj,nmasks), STAT= ierr )
148	IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate region_mask array' )
149	region_mask(:,:,:) = 0
150	region_mask = int(tmp_region_mask_real(:,:,:))
151	DEALLOCATE( tmp_region_mask_real)
152
153
154	ALLOCATE( nreg_mat(nmasks), STAT= ierr )
155	IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate nreg_mat array' )
156
157	! work out the number of regions in each mask, asssuming land is 0, and the regions are consectively numbered,
158	! without missing any number, so the number of regions is the maximum number + 1 (for land). mpp_max across the
159	! processors to get the global maxima
160	DO maskno = 1,nmasks
161	tmpint = maxval(region_mask(:,:,maskno))
162	CALL mpp_max( tmpint )
163	nreg_mat(maskno) = tmpint + 1
164	END DO
165
166	IF(lwp) THEN
167	! if writing out as binary and text, open the files.
168	IF ( ln_diaregmean_bin ) THEN
169	! Open binary for region means
170	!OPEN( UNIT=73, FILE='region_mean_timeseries.dat', FORM='UNFORMATTED', STATUS='REPLACE' )
171
172
173	CALL ctl_opn( numdct_reg_bin ,'region_mean_timeseries.dat' , 'NEW', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
174
175
176	ENDIF
177
178	IF ( ln_diaregmean_ascii ) THEN
179	! Open text files for region means
180	!OPEN( UNIT=37, FILE='region_mean_timeseries.txt', FORM='FORMATTED', STATUS='REPLACE' )
181	CALL ctl_opn( numdct_reg_txt ,'region_mean_timeseries.txt' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
182	ENDIF
183	ENDIF
184	ENDIF
185
186	END SUBROUTINE dia_regmean_init
187
188	SUBROUTINE dia_calctmb_region_mean( pinfield,pouttmb )
189	!!---------------------------------------------------------------------
190	!! * ROUTINE dia_calctmb_region_mean *
191	!!
192	!! ** Purpose : Find the Top, Mid, Bottom and Top minus Bottom fields of water Column
193	!!
194	!! ** Method :
195	!! use mbathy to find surface, mid and bottom of model levels
196	!!
197	!! History :
198	!! 3.6 ! 08-14 (E. O'Dea) Routine based on dia_wri_foam
199	!!----------------------------------------------------------------------
200	!! * Modules used
201
202	! Routine to map 3d field to top, middle, bottom
203	IMPLICIT NONE
204
205
206	! Routine arguments
207	REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN ) :: pinfield ! Input 3d field and mask
208	!REAL(wp), DIMENSION(jpi, jpj, 4 ), INTENT( OUT) :: pouttmb ! Output top, middle, bottom and surface minus bed
209	REAL(wp), DIMENSION(jpi, jpj, 3 ), INTENT( OUT) :: pouttmb ! Output top, middle, bottom and surface minus bed
210
211	! Local variables
212	INTEGER :: ji,jj,jk ! Dummy loop indices
213
214	! Local Real
215	REAL(wp) :: zmdi ! set masked values
216
217	zmdi=1.e+20 !missing data indicator for masking
218
219	! Calculate top
220	pouttmb(:,:,1) = pinfield(:,:,1)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
221
222	! Calculate middle
223	!DO jj = 1,jpj
224	! DO ji = 1,jpi
225	! jk = max(1,mbathy(ji,jj)/2)
226	! pouttmb(ji,jj,2) = pinfield(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
227	! END DO
228	!END DO
229
230	! Calculate bottom, and top minus bottom
231	DO jj = 1,jpj
232	DO ji = 1,jpi
233	jk = max(1,mbathy(ji,jj) - 1)
234	!pouttmb(ji,jj,3) = pinfield(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
235	!pouttmb(ji,jj,4) = (pouttmb(ji,jj,1) - pouttmb(ji,jj,3))tmask(ji,jj,1) + zmdi(1.0-tmask(ji,jj,1))
236	pouttmb(ji,jj,2) = pinfield(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
237	pouttmb(ji,jj,3) = (pouttmb(ji,jj,1) - pouttmb(ji,jj,2))tmask(ji,jj,1) + zmdi(1.0-tmask(ji,jj,1))
238	END DO
239	END DO
240
241	END SUBROUTINE dia_calctmb_region_mean
242
243
244	SUBROUTINE dia_regmean( kt )
245	!!----------------------------------------------------------------------
246	!! * ROUTINE dia_regmean *
247	!! ** Purpose : Produce regional mean diagnostics
248	!!
249	!! ** Method : calls dia_wri_region_mean to calculate and write the regional means for a number of variables,
250	!! (calling dia_calctmb_region_mean where necessary).
251	!!
252	!! Closes all text and binary files on last time step
253	!!
254	!!
255	!!
256	!!
257	!! History :
258	!! 3.6 ! 11-16 (J. Tinker)
259	!!
260	!!--------------------------------------------------------------------
261	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbT ! temporary T workspace
262	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbS ! temporary S workspace
263	REAL(wp) :: zmdi ! set masked values
264	INTEGER, INTENT( in ) :: kt ! ocean time-step index
265
266	REAL(wp) :: zdt ! temporary reals
267	INTEGER :: i_steps, ierr ! no of timesteps per hour, allocation error index
268	INTEGER :: maskno,jj,ji,jm,nreg ! indices of mask, i and j, and number of regions
269
270	zmdi=1.e+20 !missing data indicator for maskin
271
272	IF (ln_diaregmean) THEN
273	! If regional mean calculations required by namelist
274	! -----------------
275	! identify hourly time steps (not used)
276	zdt = rdt
277	IF( nacc == 1 ) zdt = rdtmin
278
279	IF( MOD( 3600,INT(zdt) ) == 0 ) THEN
280	i_steps = 3600/INT(zdt)
281	ELSE
282	CALL ctl_stop('STOP', 'dia_regmean: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
283	ENDIF
284
285	!i_steps = 1
286
287	!Extract 2d fields from 3d T and S with dia_calctmb_region_mean
288	CALL wrk_alloc( jpi , jpj, 3 , zwtmbT )
289	CALL wrk_alloc( jpi , jpj, 3 , zwtmbS )
290
291	CALL dia_calctmb_region_mean( tsn(:,:,:,jp_tem),zwtmbT)
292	CALL dia_calctmb_region_mean( tsn(:,:,:,jp_sal),zwtmbS)
293
294	tmp_field_mat(:,:,1) = tmp_field_mat(:,:,1) + (zwtmbT(:,:,1)*tmask(:,:,1))
295	tmp_field_mat(:,:,2) = tmp_field_mat(:,:,2) + (zwtmbT(:,:,2)*tmask(:,:,1))
296	tmp_field_mat(:,:,3) = tmp_field_mat(:,:,3) + (zwtmbT(:,:,3)*tmask(:,:,1))
297	tmp_field_mat(:,:,4) = tmp_field_mat(:,:,4) + (zwtmbS(:,:,1)*tmask(:,:,1))
298	tmp_field_mat(:,:,5) = tmp_field_mat(:,:,5) + (zwtmbS(:,:,2)*tmask(:,:,1))
299	tmp_field_mat(:,:,6) = tmp_field_mat(:,:,6) + (zwtmbS(:,:,3)*tmask(:,:,1))
300	tmp_field_mat(:,:,7) = tmp_field_mat(:,:,7) + (sshn(:,:)*tmask(:,:,1))
301
302
303	IF( ln_diaregmean_karamld ) THEN
304	tmp_field_mat(:,:,8) = tmp_field_mat(:,:,8) + (hmld_kara(:,:)*tmask(:,:,1)) !hmlp(:,:)
305	ENDIF
306	IF( ln_diaregmean_pea ) THEN
307	tmp_field_mat(:,:,9) = tmp_field_mat(:,:,9) + (pea(:,:)*tmask(:,:,1))
308	tmp_field_mat(:,:,10) = tmp_field_mat(:,:,10) + (peat(:,:)*tmask(:,:,1))
309	tmp_field_mat(:,:,11) = tmp_field_mat(:,:,11) + (peas(:,:)*tmask(:,:,1))
310	ENDIF
311
312	IF( ln_diaregmean_diaar5 ) THEN
313	tmp_field_AR5_mat(:,:,1) = tmp_field_AR5_mat(:,:,1) + (sshsteric_mat(:,:)*tmask(:,:,1))
314	tmp_field_AR5_mat(:,:,2) = tmp_field_AR5_mat(:,:,2) + (sshthster_mat(:,:)*tmask(:,:,1))
315	tmp_field_AR5_mat(:,:,3) = tmp_field_AR5_mat(:,:,3) + (sshhlster_mat(:,:)*tmask(:,:,1))
316	tmp_field_AR5_mat(:,:,4) = tmp_field_AR5_mat(:,:,4) + (zbotpres_mat(:,:)*tmask(:,:,1))
317
318	ENDIF
319
320	IF( ln_diaregmean_diasbc ) THEN
321
322	tmp_field_SBC_mat(:,:,1) = tmp_field_SBC_mat(:,:,1) + ((qsr + qns)*tmask(:,:,1))
323	tmp_field_SBC_mat(:,:,2) = tmp_field_SBC_mat(:,:,2) + (qsr*tmask(:,:,1))
324	tmp_field_SBC_mat(:,:,3) = tmp_field_SBC_mat(:,:,3) + (qns*tmask(:,:,1))
325	tmp_field_SBC_mat(:,:,4) = tmp_field_SBC_mat(:,:,4) + (emp*tmask(:,:,1))
326	tmp_field_SBC_mat(:,:,5) = tmp_field_SBC_mat(:,:,5) + (wndm*tmask(:,:,1))
327	tmp_field_SBC_mat(:,:,6) = tmp_field_SBC_mat(:,:,6) + (pressnow*tmask(:,:,1))
328	tmp_field_SBC_mat(:,:,7) = tmp_field_SBC_mat(:,:,7) + (rnf*tmask(:,:,1))
329
330
331	ENDIF
332
333	tmp_field_cnt = tmp_field_cnt + 1
334
335	IF( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 ) THEN
336
337
338	CALL dia_wri_region_mean(kt, "sst" , tmp_field_mat(:,:,1)/real(tmp_field_cnt,wp))
339	CALL dia_wri_region_mean(kt, "nbt" , tmp_field_mat(:,:,2)/real(tmp_field_cnt,wp))
340	CALL dia_wri_region_mean(kt, "dft" , tmp_field_mat(:,:,3)/real(tmp_field_cnt,wp))
341
342	CALL dia_wri_region_mean(kt, "sss" , tmp_field_mat(:,:,4)/real(tmp_field_cnt,wp))
343	CALL dia_wri_region_mean(kt, "nbs" , tmp_field_mat(:,:,5)/real(tmp_field_cnt,wp))
344	CALL dia_wri_region_mean(kt, "dfs" , tmp_field_mat(:,:,6)/real(tmp_field_cnt,wp))
345
346	CALL dia_wri_region_mean(kt, "ssh" , tmp_field_mat(:,:,7)/real(tmp_field_cnt,wp))
347
348
349	IF( ln_diaregmean_karamld ) THEN
350
351	CALL dia_wri_region_mean(kt, "mldkara" , tmp_field_mat(:,:,8)/real(tmp_field_cnt,wp)) ! tm
352	ENDIF
353	IF( ln_diaregmean_pea ) THEN
354
355	CALL dia_wri_region_mean(kt, "pea" , tmp_field_mat(:,:,9)/real(tmp_field_cnt,wp))
356	CALL dia_wri_region_mean(kt, "peat" , tmp_field_mat(:,:,10)/real(tmp_field_cnt,wp))
357	CALL dia_wri_region_mean(kt, "peas" , tmp_field_mat(:,:,11)/real(tmp_field_cnt,wp)) ! tmb
358	ENDIF
359
360	tmp_field_mat(:,:,:) = 0.
361
362	IF( ln_diaregmean_diaar5 ) THEN
363
364	CALL dia_wri_region_mean(kt, "ssh_steric" , tmp_field_AR5_mat(:,:,1)/real(tmp_field_cnt,wp))
365	CALL dia_wri_region_mean(kt, "ssh_thermosteric", tmp_field_AR5_mat(:,:,2)/real(tmp_field_cnt,wp))
366	CALL dia_wri_region_mean(kt, "ssh_halosteric" , tmp_field_AR5_mat(:,:,3)/real(tmp_field_cnt,wp))
367	CALL dia_wri_region_mean(kt, "bot_pres" , tmp_field_AR5_mat(:,:,4)/real(tmp_field_cnt,wp))
368	tmp_field_AR5_mat(:,:,:) = 0.
369	ENDIF
370
371	IF( ln_diaregmean_diasbc ) THEN
372
373	CALL dia_wri_region_mean(kt, "qt" , tmp_field_SBC_mat(:,:,1)/real(tmp_field_cnt,wp))
374	CALL dia_wri_region_mean(kt, "qsr" , tmp_field_SBC_mat(:,:,2)/real(tmp_field_cnt,wp))
375	CALL dia_wri_region_mean(kt, "qns" , tmp_field_SBC_mat(:,:,3)/real(tmp_field_cnt,wp))
376	CALL dia_wri_region_mean(kt, "emp" , tmp_field_SBC_mat(:,:,4)/real(tmp_field_cnt,wp))
377	CALL dia_wri_region_mean(kt, "wspd" , tmp_field_SBC_mat(:,:,5)/real(tmp_field_cnt,wp))
378	CALL dia_wri_region_mean(kt, "mslp" , tmp_field_SBC_mat(:,:,6)/real(tmp_field_cnt,wp))
379	CALL dia_wri_region_mean(kt, "rnf" , tmp_field_SBC_mat(:,:,7)/real(tmp_field_cnt,wp))
380	tmp_field_SBC_mat(:,:,:) = 0.
381	ENDIF
382
383	tmp_field_cnt = 0
384
385	ENDIF
386
387
388	! If on the last time step, close binary and ascii files.
389	IF( kt == nitend ) THEN
390	IF(lwp) THEN
391	IF ( ln_diaregmean_bin ) THEN
392	!Closing binary files for regional mean time series.
393	CLOSE(numdct_reg_bin)
394	ENDIF
395	IF ( ln_diaregmean_ascii ) THEN
396	!Closing text files for regional mean time series.
397	CLOSE(numdct_reg_txt)
398	ENDIF
399
400	DEALLOCATE( region_mask, nreg_mat, tmp_field_mat)
401	IF( ln_diaregmean_diaar5 ) DEALLOCATE( tmp_field_AR5_mat)
402	IF( ln_diaregmean_diasbc ) DEALLOCATE( tmp_field_SBC_mat)
403	ENDIF
404	ENDIF
405
406
407	ELSE
408	CALL ctl_warn('dia_regmean: regmean diagnostic is set to false you should not have seen this')
409	ENDIF
410
411	END SUBROUTINE dia_regmean
412
413
414	SUBROUTINE dia_wri_region_mean(kt, name, infield )
415	!!---------------------------------------------------------------------
416	!! * ROUTINE dia_tmb *
417	!!
418	!! ** Purpose : Calculate and write region mean time series for 2d arrays
419	!!
420	!! ** Method :
421	!! use
422	!!
423	!! History :
424	!! ?? ! 15/10/2015 (JTinker) Routine taken from old dia_wri_foam
425	!!----------------------------------------------------------------------
426	!! * Modules used
427	!use lib_mpp
428	!use lib_fortr
429	IMPLICIT NONE
430
431	INTEGER, INTENT(in) :: kt
432	CHARACTER (len=60) , INTENT(IN ) :: name
433	REAL(wp), DIMENSION(jpi, jpj), INTENT(IN ) :: infield ! Input 3d field and mask
434
435	! Local variables
436	INTEGER, DIMENSION(jpi, jpj) :: internal_region_mask ! Input 3d field and mask
437	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zrmet_ave,zrmet_tot,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id
438	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zrmet_out
439	REAL(wp), ALLOCATABLE, DIMENSION(:) :: ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,reg_id_mat,mask_id_mat !: region_mask
440
441	REAL(wp) :: zmdi ! set masked values
442	INTEGER :: maskno,nreg ! ocean time-step indexocean time step
443	INTEGER :: ji,jj,jk,ind,jm ! Dummy loop indices
444	INTEGER :: reg_ind_cnt ! Dummy loop indices
445
446	INTEGER :: ierr
447	REAL(wp) :: tmpreal
448	CHARACTER(LEN=180) :: FormatString,nreg_string
449	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: dummy_zrmet
450	zmdi=1.e+20 !missing data indicator for maskin
451
452	!Allocate output arrays for iomput, set to zmdi, and set a region counter = 1
453	ALLOCATE( zrmet_ave(n_regions_output), STAT= ierr )
454	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_ave array' )
455	ALLOCATE( zrmet_tot(n_regions_output), STAT= ierr )
456	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_tot array' )
457	ALLOCATE( zrmet_var(n_regions_output), STAT= ierr )
458	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_var array' )
459	ALLOCATE( zrmet_cnt(n_regions_output), STAT= ierr )
460	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_cnt array' )
461	ALLOCATE( zrmet_mask_id(n_regions_output), STAT= ierr )
462	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_mask_id array' )
463	ALLOCATE( zrmet_reg_id(n_regions_output), STAT= ierr )
464	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
465
466	ALLOCATE( zrmet_out(jpi,jpj,n_regions_output), STAT= ierr )
467	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
468
469
470	zrmet_ave(:) = zmdi
471	zrmet_tot(:) = zmdi
472	zrmet_var(:) = zmdi
473	zrmet_cnt(:) = zmdi
474	zrmet_mask_id(:) = zmdi
475	zrmet_reg_id(:) = zmdi
476
477	reg_ind_cnt = 1
478
479
480	! loop though the masks
481	DO maskno = 1,nmasks
482
483
484	! For each mask, get the number of regions (nreg), and a local copy of the region.
485	nreg = nreg_mat(maskno)
486	internal_region_mask = region_mask(:,:,maskno)
487
488	! allocate temporary stat arrays, and set to zero
489	ALLOCATE( ave_mat(nreg), STAT= ierr )
490	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ave_mat array' )
491	ALLOCATE( tot_mat(nreg), STAT= ierr )
492	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate tot_mat array' )
493	ALLOCATE( num_mat(nreg), STAT= ierr )
494	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate num_mat array' )
495	ALLOCATE( var_mat(nreg), STAT= ierr )
496	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate var_mat array' )
497	ALLOCATE( ssq_mat(nreg), STAT= ierr )
498	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ssq_mat array' )
499	ALLOCATE( cnt_mat(nreg), STAT= ierr )
500	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate cnt_mat array' )
501	ALLOCATE( reg_id_mat(nreg), STAT= ierr )
502	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate reg_id_mat array' )
503	ALLOCATE( mask_id_mat(nreg), STAT= ierr )
504	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate mask_id_mat array' )
505
506	ave_mat(:) = 0.
507	tot_mat(:) = 0.
508	num_mat(:) = 0.
509	var_mat(:) = 0.
510	cnt_mat(:) = 0.
511	ssq_mat(:) = 0.
512	reg_id_mat(:) = 0.
513	mask_id_mat(:) = 0.
514
515	! loop though the array. for each sea grid box where tmask == 1),
516	! read which region the grid box is in, add the value of the gridbox (and its square)
517	! to the total for that region, and then increment the counter for that region.
518	!CALL cpu_time(start_reg_mean_loop)
519	!WRITE(numout,*) kt,start_reg_mean_loop
520	DO ji = nldi,nlei
521	DO jj = nldj,nlej
522	IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
523	ind = internal_region_mask(ji,jj)+1
524	tot_mat(ind) = tot_mat(ind) + (infield(ji,jj))
525	ssq_mat(ind) = ssq_mat(ind) + ( infield(ji,jj) * infield(ji,jj))
526	cnt_mat(ind) = cnt_mat(ind) + 1.
527	ENDIF
528	END DO
529	END DO
530	! sum the totals, the counts, and the squares across the processors
531	CALL mpp_sum( tot_mat,nreg )
532	CALL mpp_sum( ssq_mat,nreg )
533	CALL mpp_sum( cnt_mat,nreg )
534
535
536	!calculate the mean and variance from the total, sum of squares and the count.
537
538	ave_mat = tot_mat(:)/cnt_mat(:)
539	var_mat = ssq_mat(:)/cnt_mat(:) - (ave_mat(:)*ave_mat(:))
540
541
542	!mask array of mask and region number.
543	DO jj = 1,nreg
544	reg_id_mat(jj) = real(jj-1)
545	mask_id_mat(jj) = real(maskno)
546	END DO
547
548
549	!write text and binary, and note region statistics for current mask for later iom_put
550	IF( lwp ) THEN
551
552	!Write out ascii and binary if requred
553	IF ( ln_diaregmean_bin ) THEN
554	!Writing out regional mean time series to binary files
555	WRITE(numdct_reg_bin) name,kt,maskno,n_regions_output
556	WRITE(numdct_reg_bin) ave_mat
557	WRITE(numdct_reg_bin) tot_mat
558	WRITE(numdct_reg_bin) var_mat
559	WRITE(numdct_reg_bin) ssq_mat
560	WRITE(numdct_reg_bin) cnt_mat
561	ENDIF
562
563	IF ( ln_diaregmean_ascii ) THEN
564	!Writing out regional mean time series to text files
565
566	WRITE(nreg_string, "(I5)") nreg
567	FormatString = "(A17,"//trim(nreg_string)//"F15.3)"
568	WRITE(numdct_reg_txt, FMT="(A17,I6,I6)") name,kt,maskno
569	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(name)//" "//"ave_mat:", ave_mat
570	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(name)//" "//"tot_mat:", tot_mat
571	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(name)//" "//"var_mat:", var_mat
572	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(name)//" "//"ssq_mat:", ssq_mat
573	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(name)//" "//"cnt_mat:", cnt_mat
574	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(name)//" "//"reg_mat:", reg_id_mat
575	WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(name)//" "//"msk_mat:", mask_id_mat
576
577	ENDIF
578
579	DO jm = 1,nreg
580	zrmet_ave( reg_ind_cnt) = ave_mat(jm)
581	zrmet_tot( reg_ind_cnt) = tot_mat(jm)
582	zrmet_var( reg_ind_cnt) = var_mat(jm)
583	zrmet_cnt( reg_ind_cnt) = cnt_mat(jm)
584	zrmet_reg_id( reg_ind_cnt) = reg_id_mat(jm)
585	zrmet_mask_id(reg_ind_cnt) = mask_id_mat(jm)
586
587	reg_ind_cnt = reg_ind_cnt + 1
588	END DO
589
590	ENDIF
591
592	DEALLOCATE(ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,reg_id_mat,mask_id_mat)
593	!DEALLOCATE(tot_mat_2)
594	END DO
595
596
597	!With current field_def.xml and iodef.xml, these fields must be output, so set to dummy values if not required.
598
599	IF ( ln_diaregmean_nc ) THEN
600
601
602
603	DO jm = 1,n_regions_output
604	zrmet_out(:,:,jm) = zrmet_ave(jm)
605	END DO
606	CALL iom_put( "reg_" // trim(name) // '_ave', zrmet_out )
607	zrmet_out(:,:,:) = 0
608
609	DO jm = 1,n_regions_output
610	zrmet_out(:,:,jm) = zrmet_tot(jm)
611	END DO
612	CALL iom_put( "reg_" // trim(name) // '_tot', zrmet_out )
613	zrmet_out(:,:,:) = 0
614
615	DO jm = 1,n_regions_output
616	zrmet_out(:,:,jm) = zrmet_var(jm)
617	END DO
618	CALL iom_put( "reg_" // trim(name) // '_var', zrmet_out )
619	zrmet_out(:,:,:) = 0
620
621	DO jm = 1,n_regions_output
622	zrmet_out(:,:,jm) = zrmet_cnt(jm)
623	END DO
624	CALL iom_put( "reg_" // trim(name) // '_cnt', zrmet_out )
625	zrmet_out(:,:,:) = 0
626
627	DO jm = 1,n_regions_output
628	zrmet_out(:,:,jm) = zrmet_reg_id(jm)
629	END DO
630	CALL iom_put( "reg_" // trim(name) // '_reg_id', zrmet_out )
631	zrmet_out(:,:,:) = 0
632
633	DO jm = 1,n_regions_output
634	zrmet_out(:,:,jm) = zrmet_mask_id(jm)
635	END DO
636	CALL iom_put( "reg_" // trim(name) // '_mask_id', zrmet_out )
637	zrmet_out(:,:,:) = 0
638	ELSE
639
640	ALLOCATE( dummy_zrmet(jpi,jpj,n_regions_output), STAT= ierr )
641	IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate dummy_zrmet array' )
642
643	DO jm = 1,n_regions_output
644	dummy_zrmet(:,:,jm) = real(jm,wp)
645	END DO
646
647	DO jm = 1,9
648	CALL iom_put( "reg_" // trim(name) // '_ave', dummy_zrmet )
649	CALL iom_put( "reg_" // trim(name) // '_tot', dummy_zrmet )
650	CALL iom_put( "reg_" // trim(name) // '_var', dummy_zrmet )
651	CALL iom_put( "reg_" // trim(name) // '_cnt', dummy_zrmet )
652	CALL iom_put( "reg_" // trim(name) // '_reg_id', dummy_zrmet )
653	CALL iom_put( "reg_" // trim(name) // '_mask_id', dummy_zrmet )
654	END DO
655
656	DEALLOCATE( dummy_zrmet)
657	ENDIF
658
659	DEALLOCATE(zrmet_ave,zrmet_tot,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id,zrmet_out)
660	END SUBROUTINE dia_wri_region_mean
661
662
663	!!======================================================================
664	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/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90 @ 10489

Download in other formats: