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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diaregmean.F90 in branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90 @ 7613

Last change on this file since 7613 was 7613, checked in by hadjt, 7 years ago
Added namelist entries to effectively switch off region mean, pea and NOOS code. Also use correct nemo file open subroutines, rahter than low level open statements with fixed file ids (37 and 73)
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 = 1,jpi
521	DO jj = 1,jpj
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.

Download in other formats: