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

Last change on this file since 10960 was 10960, checked in by rrenshaw, 19 months ago

code fix

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