Context Navigation

diaregmean.F90 @ 11038

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

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

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

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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90 @ 11038

Download in other formats: