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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diaregmean.F90 in NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diaregmean.F90 @ 15397

Last change on this file since 15397 was 15397, checked in by hadjt, 12 months ago

Region mean and NOOS transport XIOS scaler output
Scalar output, so passing 85 and 3 values per iom call rather than 297*375*85.
DIADCT edited (nd changes commmeted out) to allow compilation on IFORT.
In all NEMO, reading section_ijglobal.diadct was the only compilation error, ifort didn't like ctl_opn, and cray didn't like call ftell.

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