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 @ 15411

Last change on this file since 15411 was 15411, checked in by hadjt, 11 months ago

Region means updated to calculate river and precip/evap heat flux given surface temperature

File size: 67.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,9),  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(9)  ::    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(9)             ::    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
694            tmp_field_SBC_mat(:,:,8) = tmp_field_SBC_mat(:,:,8) + (emp*tmask(:,:,1)*tsn(:,:,1,jp_tem)*3850.)
695            name_SBC_mat(8) = 'empheat'
696            tmp_field_SBC_mat(:,:,9) = tmp_field_SBC_mat(:,:,9) + (rnf*tmask(:,:,1)*tsn(:,:,1,jp_tem)*3850.)
697            name_SBC_mat(9) = 'rnfheat'
698
699        ENDIF
700
701        output_mulitpler_dat_mat(:) = 1.
702        output_mulitpler_dat_mat(6)  = output_mulitpler_HSVM_mat(1) ! 1e-12
703        output_mulitpler_dat_mat(12) = output_mulitpler_HSVM_mat(2) ! 1e-12
704        output_mulitpler_AR5_mat(:) = 1.
705        output_mulitpler_SBC_mat(:) = 1.
706
707        ! On the hour, calculate hourly means from the hourly total,and process the regional means.
708
709        tmp_field_cnt = tmp_field_cnt + 1
710
711       
712        IF ( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 ) THEN
713
714           
715            DO vi=1,19 ! State loop
716
717               do_reg_mean = .TRUE.
718
719               IF (vi == 16) THEN
720                 IF( .not. ln_diaregmean_mld ) do_reg_mean = .FALSE.   
721               ENDIF
722
723               IF ((vi == 17) .OR. (vi == 18) .OR. (vi == 19) ) THEN
724                 IF( .not. ln_diaregmean_pea ) do_reg_mean = .FALSE.   
725               ENDIF
726
727               tmp_name=TRIM( name_dat_mat(vi) )
728               IF ( do_reg_mean ) THEN
729                   IF (iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_ave'))))    .OR. &
730                     & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_tot'))))    .OR. &
731                     & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_totheat'))))    .OR. &
732                     & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_var'))))    .OR. &
733                     & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))))    .OR. &
734                     & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_reg_id')))) .OR. &
735                     & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_mask_id')))) ) THEN
736
737                        CALL dia_wri_region_mean(kt, TRIM(tmp_name) , output_mulitpler_dat_mat(vi)*tmp_field_mat(:,:,vi)/real(tmp_field_cnt,wp))
738                        !WRITE(numout,*)  'JT dia_regmean SBC variable - region mean: ',TRIM( name_dat_mat(vi) ),';'
739                    !ELSE
740                       !WRITE(numout,*)  'JT dia_regmean SBC variable - no iom_use: ',TRIM( name_dat_mat(vi) ),';'
741                    ENDIF
742                !ELSE
743                    !WRITE(numout,*)  'JT dia_regmean SBC variable - no do_reg_mean: ',TRIM( name_dat_mat(vi) ),';',ln_diaregmean_mld,ln_diaregmean_pea
744                ENDIF
745                tmp_name=""
746            END DO
747           
748            tmp_field_mat(:,:,:) = 0.
749
750            DO vi=1,4 ! State loop
751
752                tmp_name=TRIM( name_HSCM_mat(vi) ) // trim('_inst')
753                IF (iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_ave'))))    .OR. &
754                  & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_tot'))))    .OR. &
755                  & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_totarea'))))    .OR. &
756                  & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_var'))))    .OR. &
757                  & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))))    .OR. &
758                  & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_reg_id')))) .OR. &
759                  & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_mask_id')))) ) THEN
760
761                    CALL dia_wri_region_mean(kt, TRIM(tmp_name) , output_mulitpler_HSVM_mat(vi)*tmp_field_HSVM_mat(:,:,vi))
762                ENDIF
763                tmp_name=""
764            END DO
765
766            tmp_field_HSVM_mat(:,:,:) = 0.
767            IF( ln_diaregmean_diaar5  ) THEN
768                DO vi=1,4 ! State loop
769
770                    tmp_name=TRIM( name_AR5_mat(vi) )
771                    IF (iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_ave'))))    .OR. &
772                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_tot'))))    .OR. &
773                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_totarea'))))    .OR. &
774                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_var'))))    .OR. &
775                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))))    .OR. &
776                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_reg_id')))) .OR. &
777                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_mask_id')))) ) THEN
778
779                        CALL dia_wri_region_mean(kt, TRIM(tmp_name) , output_mulitpler_AR5_mat(vi)*tmp_field_AR5_mat(:,:,vi)/real(tmp_field_cnt,wp))
780                    ENDIF
781                    tmp_name=""
782                END DO
783                tmp_field_AR5_mat(:,:,:) = 0.
784            ENDIF
785
786            IF( ln_diaregmean_diasbc  ) THEN
787                DO vi=1,9 ! State loop
788
789                    tmp_name=TRIM( name_SBC_mat(vi) )
790                    IF (iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_ave'))))    .OR. &
791                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_tot'))))    .OR. &
792                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_totarea'))))    .OR. &
793                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_var'))))    .OR. &
794                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))))    .OR. &
795                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_reg_id')))) .OR. &
796                      & iom_use(trim( trim(trim("reg_") // trim(tmp_name) // trim('_mask_id')))) ) THEN
797
798                        CALL dia_wri_region_mean(kt, TRIM(tmp_name) , output_mulitpler_SBC_mat(vi)*tmp_field_SBC_mat(:,:,vi)/real(tmp_field_cnt,wp))
799                    ENDIF
800                    tmp_name=""
801                END DO
802                tmp_field_SBC_mat(:,:,:) = 0.
803            ENDIF
804
805
806#if defined key_fabm
807            !ADD Biogeochemistry
808           
809            IF( ln_diaregmean_bgc  ) THEN !ln_diaregmean_bgc
810               
811                ! Loop through 3d BGC tracers
812                DO jn=1,jp_fabm ! State loop
813
814                    ! get variable name for different levels
815                    tmp_name_bgc_top=TRIM(TRIM("top_")//TRIM(ctrcnm(jn)))
816                    tmp_name_bgc_bot=TRIM(TRIM("bot_")//TRIM(ctrcnm(jn)))
817                    tmp_name_bgc_dif=TRIM(TRIM("dif_")//TRIM(ctrcnm(jn)))
818                    tmp_name_bgc_zav=TRIM(TRIM("zav_")//TRIM(ctrcnm(jn)))
819                    tmp_name_bgc_vol=TRIM(TRIM("vol_")//TRIM(ctrcnm(jn)))
820
821                    ! print out names if verbose
822                    IF(verbose .AND. lwp) THEN
823                        WRITE(numout,*)
824                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_top : ',TRIM(tmp_name_bgc_top)
825                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_bot : ',TRIM(tmp_name_bgc_bot)
826                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_dif : ',TRIM(tmp_name_bgc_dif)
827                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_zav : ',TRIM(tmp_name_bgc_zav)
828                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_vol : ',TRIM(tmp_name_bgc_vol)
829                        CALL FLUSH(numout) 
830
831                    ENDIF
832
833                    !Allocate working array, and get surface, bed etc fields.
834                    !CALL wrk_alloc( jpi , jpj,  6 , zwtmbBGC )
835                    ALLOCATE (zwtmbBGC(jpi , jpj, 6),  STAT= ierr )
836                    IF( ierr /= 0 )   CALL ctl_stop( 'dia_regmean: failed to allocate zwtmbBGC array' )
837                    CALL dia_calctmb_region_mean(  trn(:,:,:,jn),zwtmbBGC )
838
839
840                    !Print out 2d fields to ascii text files to check values if verbose. (24MB per time step, per BGC variable)
841                    IF (verbose) THEN
842
843                        WRITE (tmp_output_filename, "(A4,I3.3,A1,I6.6,A1,I3.3,A4)") "bgc_",jn,"_",kt,"_",narea,".txt"
844                        WRITE (*,*) tmp_output_filename
845                        OPEN(UNIT=74,FILE=TRIM(tmp_output_filename))
846
847                        DO ji = 1,jpi
848                            DO jj = 1,jpj
849                                WRITE(74,FMT="(I4,I4,F3,F25.5,F25.5,F25.5,F25.5,F25.5)") nimpp+ji, njmpp+jj,tmask(ji,jj,1),&
850                                      & zwtmbBGC(ji,jj,1),zwtmbBGC(ji,jj,2),zwtmbBGC(ji,jj,3),zwtmbBGC(ji,jj,4),zwtmbBGC(ji,jj,5)/1e9
851                            END DO
852                        END DO
853                        CLOSE(74)
854                    ENDIF
855                   
856                    ! Do region means
857                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_top)   , zwtmbBGC(:,:,1))
858                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_bot)   , zwtmbBGC(:,:,2))
859                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_dif)   , zwtmbBGC(:,:,3))
860                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_zav)   , zwtmbBGC(:,:,4))
861                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_vol)   , zwtmbBGC(:,:,5)/1e9)
862
863
864                    !Deallocate working array
865                    !JT CALL wrk_dealloc( jpi , jpj,  6 , zwtmbBGC )
866                    DEALLOCATE ( zwtmbBGC )
867                ENDDO ! State loop
868            ENDIF !ln_diaregmean_bgc
869
870#endif
871             
872            tmp_field_cnt = 0
873 
874        ENDIF ! ( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 )
875       
876       
877        ! If on the last time step, close binary and ascii files.
878        IF( kt == nitend ) THEN
879            IF(lwp) THEN
880                IF ( ln_diaregmean_bin ) THEN
881                    !Closing binary files for regional mean time series.
882                    CLOSE(numdct_reg_bin)
883                ENDIF
884                IF ( ln_diaregmean_ascii ) THEN
885                    !Closing text files for regional mean time series.
886                    CLOSE(numdct_reg_txt)
887                ENDIF
888
889                DEALLOCATE( region_mask, nreg_mat, tmp_field_mat,tmp_field_HSVM_mat)
890                IF( ln_diaregmean_diaar5  ) DEALLOCATE( tmp_field_AR5_mat)
891                IF( ln_diaregmean_diasbc  ) DEALLOCATE( tmp_field_SBC_mat)
892            ENDIF
893        ENDIF
894         
895         
896      ELSE
897        CALL ctl_warn('dia_regmean: regmean diagnostic is set to false you should not have seen this')
898      ENDIF
899     
900   END SUBROUTINE dia_regmean
901   
902   
903   SUBROUTINE dia_wri_region_mean(kt, tmp_name,         infield  )
904      !!---------------------------------------------------------------------
905      !!                  ***  ROUTINE dia_tmb  ***
906      !!                   
907      !! ** Purpose :   Calculate and write region mean time series for 2d arrays
908      !!
909      !! ** Method  :   
910      !!      use
911      !!
912      !! History :
913      !!   ??  !  15/10/2015  (JTinker) Routine taken from old dia_wri_foam
914      !!----------------------------------------------------------------------
915      !! * Modules used
916      !use lib_mpp
917      !use lib_fortr
918      IMPLICIT NONE
919     
920      INTEGER, INTENT(in) ::   kt
921      CHARACTER (len=*) , INTENT(IN   ) ::    tmp_name
922      REAL(wp), DIMENSION(jpi, jpj), INTENT(IN   ) :: infield    ! Input 3d field and mask
923     
924      ! Local variables
925      INTEGER, DIMENSION(jpi, jpj) :: internal_region_mask    ! Input 3d field and mask
926      REAL(wp), DIMENSION(jpi, jpj) :: internal_infield    ! Internal data field
927      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zrmet_ave,zrmet_tot,zrmet_totarea,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id  ,zrmet_min,zrmet_max
928      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zrmet_out
929      REAL(wp), ALLOCATABLE,   DIMENSION(:) ::   ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,reg_id_mat,mask_id_mat,area_mat,totarea_mat    !: region_mask
930      !REAL(wp), ALLOCATABLE,   DIMENSION(:) ::   min_mat,max_mat   !: region_mask
931     
932      REAL(wp)                         ::   zmdi, zrmet_val      ! set masked values
933      INTEGER :: maskno,nreg  ! ocean time-step indexocean time step           
934      INTEGER :: ji,jj,jk,ind,jm ! Dummy loop indices
935      INTEGER :: reg_ind_cnt ! Dummy loop indices
936     
937      INTEGER  ::   ierr     
938      REAL(wp)  :: tmpreal
939      CHARACTER(LEN=180) :: FormatString,nreg_string,tmp_name_iom
940      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   dummy_zrmet
941      LOGICAL       ::   verbose     
942
943
944      verbose = ln_diaregmean_verbose
945
946
947      zmdi=1.e+20 !missing data indicator for maskin
948     
949      !Allocate output arrays for iomput, set to zmdi, and set a region counter = 1
950      ALLOCATE( zrmet_ave(n_regions_output),  STAT= ierr )
951        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_ave array' )
952      ALLOCATE( zrmet_tot(n_regions_output),  STAT= ierr )
953        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_tot array' )
954      ALLOCATE( zrmet_totarea(n_regions_output),  STAT= ierr )
955        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_totarea array' )
956      ALLOCATE( zrmet_var(n_regions_output),  STAT= ierr )
957        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_var array' )
958      ALLOCATE( zrmet_cnt(n_regions_output),  STAT= ierr )
959        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_cnt array' )
960      ALLOCATE( zrmet_mask_id(n_regions_output),  STAT= ierr )
961        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_mask_id array' )
962      ALLOCATE( zrmet_reg_id(n_regions_output),  STAT= ierr )
963        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
964
965
966      ALLOCATE( zrmet_min(n_regions_output),  STAT= ierr )
967        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_min array' )
968      ALLOCATE( zrmet_max(n_regions_output),  STAT= ierr )
969        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_max array' )
970     
971      ALLOCATE( zrmet_out(n_regions_output),  STAT= ierr )
972        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
973
974 
975     
976        IF(lwp .AND. verbose) THEN
977              WRITE(numout,*)
978              WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//';'
979              WRITE(numout,*)
980        ENDIF
981       
982      DO ji = 1,jpi
983          DO jj = 1,jpj
984            internal_infield(ji,jj) = infield(ji,jj)
985          END DO
986      END DO
987       
988      ! Check for NANS # JT   03/09/2018
989      DO ji = 1,jpi
990          DO jj = 1,jpj
991              IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
992                  IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
993                      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)
994                      internal_infield(ji,jj) = 0.
995                  ENDIF
996              ELSE               
997                  IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
998                      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)
999                      internal_infield(ji,jj) = 0.
1000                  ENDIF
1001              ENDIF
1002          END DO
1003      END DO
1004     
1005     
1006      zrmet_ave(:) = zmdi
1007      zrmet_tot(:) = zmdi
1008      zrmet_totarea(:) = zmdi
1009      zrmet_var(:) = zmdi
1010      zrmet_cnt(:) = zmdi
1011      zrmet_mask_id(:) = zmdi
1012      zrmet_reg_id(:) = zmdi
1013     
1014      zrmet_min(:) = zmdi
1015      zrmet_max(:) = zmdi
1016      reg_ind_cnt = 1
1017     
1018     
1019      ! loop though the masks
1020      DO maskno = 1,nmasks
1021          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin mask loops: ',maskno
1022         
1023         
1024          ! For each mask, get the number of regions (nreg), and a local copy of the region.
1025          nreg = nreg_mat(maskno)
1026          internal_region_mask = region_mask(:,:,maskno)
1027         
1028          ! allocate temporary stat arrays, and set to zero
1029          ALLOCATE( ave_mat(nreg),  STAT= ierr )
1030          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ave_mat array' )
1031          ALLOCATE( tot_mat(nreg),      STAT= ierr )
1032          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate tot_mat array' )
1033          ALLOCATE( num_mat(nreg),  STAT= ierr )
1034          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate num_mat array' )
1035          ALLOCATE( var_mat(nreg),  STAT= ierr )
1036          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate var_mat array' )
1037          ALLOCATE( ssq_mat(nreg),  STAT= ierr )
1038          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ssq_mat array' )
1039          ALLOCATE( cnt_mat(nreg),  STAT= ierr )
1040          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate cnt_mat array' )
1041          ALLOCATE( area_mat(nreg),  STAT= ierr )
1042          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate area_mat array' )
1043          ALLOCATE( totarea_mat(nreg),  STAT= ierr )
1044          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate totarea_mat array' )
1045
1046          !ALLOCATE( min_mat(nreg),  STAT= ierr )
1047          !IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate min_mat array' )
1048          !ALLOCATE( max_mat(nreg),  STAT= ierr )
1049          !IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate max_mat array' )
1050
1051          ALLOCATE( reg_id_mat(nreg),  STAT= ierr )
1052          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate reg_id_mat array' )
1053          ALLOCATE( mask_id_mat(nreg),  STAT= ierr )
1054          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate mask_id_mat array' )
1055
1056
1057         
1058          ave_mat(:) = 0.
1059          tot_mat(:) = 0.
1060          num_mat(:) = 0.
1061          var_mat(:) = 0.
1062          cnt_mat(:) = 0.
1063          ssq_mat(:) = 0.
1064          area_mat(:) = 0.
1065          totarea_mat(:) = 0.
1066
1067
1068          !min_mat(:) = zmdi
1069          !max_mat(:) = -zmdi
1070          reg_id_mat(:) = 0.
1071          mask_id_mat(:) = 0.
1072         
1073          ! loop though the array. for each sea grid box where tmask == 1),
1074          ! read which region the grid box is in, add the value of the gridbox (and its square)
1075          ! to the total for that region, and then increment the counter for that region.
1076          !CALL cpu_time(start_reg_mean_loop)
1077          !WRITE(numout,*) kt,start_reg_mean_loop
1078          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin spatial loops: '
1079          DO ji = nldi,nlei
1080              DO jj = nldj,nlej
1081                    IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
1082                        ind = internal_region_mask(ji,jj)+1
1083                        !tot_mat(ind) = tot_mat(ind) + (internal_infield(ji,jj))
1084                        !ssq_mat(ind) = ssq_mat(ind) + ( internal_infield(ji,jj) *  internal_infield(ji,jj))
1085                        !cnt_mat(ind) = cnt_mat(ind) + 1.
1086                        ! Area Weighted values - region_area_mat == 1. or area depending on ln_diaregmean_areawgt
1087                        totarea_mat(ind) = totarea_mat(ind) + (region_area_mat(ji,jj)*internal_infield(ji,jj))
1088                        tot_mat(ind) = tot_mat(ind) + (internal_infield(ji,jj))
1089                        ssq_mat(ind) = ssq_mat(ind) + (region_area_mat(ji,jj)*(internal_infield(ji,jj) * internal_infield(ji,jj)))
1090                        cnt_mat(ind) = cnt_mat(ind) + 1.
1091                        area_mat(ind) = area_mat(ind) + (region_area_mat(ji,jj)*1.)
1092
1093
1094
1095                        !min_mat(ind) = min(min_mat(ind),internal_infield(ji,jj))
1096                        !max_mat(ind) = max(max_mat(ind),internal_infield(ji,jj))
1097                    ENDIF
1098              END DO
1099          END DO
1100          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finish spatial loops: '
1101          ! sum the totals, the counts, and the squares across the processors         
1102          CALL mpp_sum( 'diaregionmean',tot_mat,nreg )
1103          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum tot'
1104          CALL mpp_sum( 'diaregionmean',cnt_mat,nreg )
1105          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum cnt'
1106          CALL mpp_sum( 'diaregionmean',area_mat,nreg )
1107          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum area'
1108          CALL mpp_sum( 'diaregionmean',totarea_mat,nreg )
1109          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum totarea_mat'
1110
1111
1112
1113          !tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_var'))
1114          !IF (iom_use(trim(tmp_name_iom)) .OR. ln_diaregmean_bin .OR. ln_diaregmean_ascii) THEN
1115              CALL mpp_sum( 'diaregionmean',ssq_mat,nreg )
1116              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum ssq'
1117          !ENDIF
1118         
1119   
1120
1121          !tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_min'))
1122          !IF (iom_use(trim(tmp_name_iom)) .OR. ln_diaregmean_bin .OR. ln_diaregmean_ascii) THEN
1123              !CALL mpp_min( 'diaregionmean',min_mat,nreg )
1124              !IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_min'
1125          !ENDIF
1126         
1127
1128          !tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_max'))
1129          !IF (iom_use(trim(tmp_name_iom)) .OR. ln_diaregmean_bin .OR. ln_diaregmean_ascii)  THEN
1130              !CALL mpp_max( 'diaregionmean',max_mat,nreg )
1131              !IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_max'
1132          !ENDIF
1133         
1134          ! calculate the mean and variance from the total, sum of squares and the count.
1135          ! When area weighting, you can't area weight the total.
1136          ! this if block may be redundant, as totarea_mat == tot_mat, and cnt_mat == area_mat when ln_diaregmean_areawgt == False
1137          IF (ln_diaregmean_areawgt) THEN
1138            ave_mat = totarea_mat(:)/area_mat(:)
1139            var_mat = ssq_mat(:)/area_mat(:) - (ave_mat(:)*ave_mat(:))
1140          ELSE
1141            ave_mat = tot_mat(:)/cnt_mat(:)
1142            var_mat = ssq_mat(:)/cnt_mat(:) - (ave_mat(:)*ave_mat(:))
1143          ENDIF
1144         
1145         
1146         
1147          !mask array of mask and region number.
1148          DO jj = 1,nreg
1149              reg_id_mat(jj) = real(jj-1)
1150              mask_id_mat(jj) = real(maskno)
1151          END DO
1152         
1153         
1154          !write text and binary, and note region statistics for current mask for later iom_put
1155          IF( lwp ) THEN 
1156         
1157              !Write out ascii and binary if requred
1158              IF ( ln_diaregmean_bin ) THEN
1159                  !Writing out regional mean time series to binary files
1160                  WRITE(numdct_reg_bin) tmp_name,kt,maskno,n_regions_output
1161                  WRITE(numdct_reg_bin) ave_mat
1162                  WRITE(numdct_reg_bin) tot_mat
1163                  WRITE(numdct_reg_bin) var_mat
1164                  WRITE(numdct_reg_bin) ssq_mat
1165                  WRITE(numdct_reg_bin) cnt_mat
1166                  !WRITE(numdct_reg_bin) min_mat
1167                  !WRITE(numdct_reg_bin) max_mat
1168              ENDIF
1169             
1170              IF ( ln_diaregmean_ascii  ) THEN
1171                  !Writing out regional mean time series to text files
1172
1173                  WRITE(nreg_string, "(I5)") nreg
1174                  FormatString = "(A30,"//trim(nreg_string)//"F25.3)"
1175                  WRITE(numdct_reg_txt, FMT="(A30,I6,I6)") tmp_name,kt,maskno           
1176                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ave_mat:", ave_mat
1177                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"tot_mat:", tot_mat
1178                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"var_mat:", var_mat
1179                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ssq_mat:", ssq_mat
1180                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"cnt_mat:", cnt_mat
1181                  !WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"min_mat:", min_mat
1182                  !WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"max_mat:", max_mat
1183                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"reg_mat:", reg_id_mat
1184                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"msk_mat:", mask_id_mat
1185
1186              ENDIF
1187          ENDIF
1188             
1189          ! JT Fixed, was not meant to be inside the lwp if block
1190          DO jm = 1,nreg
1191              zrmet_ave(    reg_ind_cnt) =     ave_mat(jm)
1192              zrmet_tot(    reg_ind_cnt) =     tot_mat(jm)
1193              zrmet_totarea(    reg_ind_cnt) =     totarea_mat(jm)
1194              zrmet_var(    reg_ind_cnt) =     var_mat(jm)
1195              zrmet_cnt(    reg_ind_cnt) =     cnt_mat(jm)
1196              !zrmet_min(    reg_ind_cnt) =     min_mat(jm)
1197              !zrmet_max(    reg_ind_cnt) =     max_mat(jm)
1198              zrmet_reg_id( reg_ind_cnt) =  reg_id_mat(jm)
1199              zrmet_mask_id(reg_ind_cnt) = mask_id_mat(jm)
1200           
1201              reg_ind_cnt = reg_ind_cnt + 1 
1202          END DO
1203     
1204       
1205          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean about to deallocated arrays for ',kt,maskno
1206          DEALLOCATE(ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,reg_id_mat,mask_id_mat,totarea_mat, area_mat)
1207          !DEALLOCATE(min_mat,max_mat)
1208
1209          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean deallocated arrays for ',kt,maskno
1210          IF(lwp .AND. ln_diaregmean_ascii ) CALL FLUSH(numdct_reg_txt)
1211          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean flushed region mean text for ',kt,maskno
1212      END DO
1213
1214      IF(lwp .AND. verbose) THEN                   ! Control print
1215         WRITE(numout,*) 'dia_regmean ready to start iom_put'
1216         CALL FLUSH(numout)
1217      ENDIF
1218     
1219      !With current field_def.xml and iodef.xml, these fields must be output, so set to dummy values if not required.
1220     
1221      IF ( ln_diaregmean_nc ) THEN
1222     
1223          zrmet_out(:) = 0
1224          zrmet_val = 0
1225          tmp_name_iom = ''
1226
1227          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean ready to start iom_put: ',trim(tmp_name)
1228         
1229          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_ave'))
1230          IF (iom_use(trim(tmp_name_iom))) THEN
1231              DO jm = 1,n_regions_output
1232                zrmet_val = zrmet_ave(jm)
1233    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1234    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1235                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1236                zrmet_out(jm) = zrmet_val
1237              END DO     
1238              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1239              CALL iom_put(trim(tmp_name_iom), zrmet_out(:) ) 
1240              zrmet_out(:) = 0
1241              zrmet_val = 0
1242              tmp_name_iom = ''
1243          ENDIF
1244
1245          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_tot'))
1246          IF (iom_use(trim(tmp_name_iom))) THEN
1247              DO jm = 1,n_regions_output
1248                zrmet_val = zrmet_tot(jm)
1249    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1250    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1251                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1252                zrmet_out(jm) = zrmet_val
1253              END DO
1254              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1255              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) ) 
1256              zrmet_out(:) = 0
1257              zrmet_val = 0
1258              tmp_name_iom = ''
1259          ENDIF
1260
1261
1262
1263          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_totarea'))
1264          IF (iom_use(trim(tmp_name_iom))) THEN
1265              DO jm = 1,n_regions_output
1266                zrmet_val = zrmet_totarea(jm)
1267    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1268    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1269                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1270                zrmet_out(jm) = zrmet_val
1271              END DO
1272              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1273              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) ) 
1274              zrmet_out(:) = 0
1275              zrmet_val = 0
1276              tmp_name_iom = ''
1277          ENDIF
1278
1279          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_var'))
1280          IF (iom_use(trim(tmp_name_iom))) THEN
1281              DO jm = 1,n_regions_output
1282                zrmet_val = zrmet_var(jm)
1283    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1284    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1285                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1286                zrmet_out(jm) = zrmet_val
1287              END DO
1288              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1289              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) )
1290              zrmet_out(:) = 0
1291              zrmet_val = 0
1292              tmp_name_iom = ''
1293          ENDIF
1294
1295          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))
1296          IF (iom_use(trim(tmp_name_iom))) THEN
1297              DO jm = 1,n_regions_output
1298                zrmet_val = zrmet_cnt(jm)
1299    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1300    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1301                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1302                zrmet_out(jm) = zrmet_val
1303              END DO
1304              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1305              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) )
1306              zrmet_out(:) = 0
1307              zrmet_val = 0
1308              tmp_name_iom = ''
1309          ENDIF
1310
1311          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_reg_id'))
1312          IF (iom_use(trim(tmp_name_iom))) THEN
1313              DO jm = 1,n_regions_output
1314                zrmet_val = zrmet_reg_id(jm)
1315    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1316    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1317                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1318                zrmet_out(jm) = zrmet_val
1319              END DO
1320              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1321              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) ) 
1322              zrmet_out(:) = 0
1323              zrmet_val = 0
1324              tmp_name_iom = ''
1325          ENDIF
1326
1327          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_mask_id'))
1328          IF (iom_use(trim(tmp_name_iom))) THEN
1329              DO jm = 1,n_regions_output
1330                zrmet_val = zrmet_mask_id(jm)
1331    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1332    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1333                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1334                zrmet_out(jm) = zrmet_val
1335              END DO
1336              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)         
1337              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) )
1338              zrmet_out(:) = 0
1339              zrmet_val = 0
1340              tmp_name_iom = ''
1341          ENDIF
1342     
1343!      ELSE
1344!       
1345!          ALLOCATE( dummy_zrmet(jpi,jpj,n_regions_output),  STAT= ierr )
1346!            IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate dummy_zrmet array' )
1347
1348!          DO jm = 1,n_regions_output
1349!              dummy_zrmet(:,:,jm) =     real(jm,wp)
1350!          END DO
1351
1352!          DO jm = 1,9
1353!              !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 )
1354!              !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 )
1355!              !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 )
1356!              !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 )
1357!              !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 )
1358!              !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 )
1359
1360
1361!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_ave'))
1362!              IF (iom_use(trim(tmp_name_iom))) THEN
1363!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1364!              ENDIF
1365!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_tot'))
1366!              IF (iom_use(trim(tmp_name_iom))) THEN
1367!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1368!              ENDIF
1369!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_var'))
1370!              IF (iom_use(trim(tmp_name_iom))) THEN
1371!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1372!              ENDIF
1373!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))
1374!              IF (iom_use(trim(tmp_name_iom))) THEN
1375!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1376!              ENDIF
1377!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_reg_id'))
1378!              IF (iom_use(trim(tmp_name_iom))) THEN
1379!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1380!              ENDIF
1381!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_mask_id'))
1382!              IF (iom_use(trim(tmp_name_iom))) THEN
1383!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1384!              ENDIF
1385!
1386!          END DO
1387!   
1388!          DEALLOCATE( dummy_zrmet)
1389      ENDIF
1390     
1391      DEALLOCATE(zrmet_ave,zrmet_tot,zrmet_totarea,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id,zrmet_min,zrmet_max,zrmet_out)
1392
1393      IF(lwp .AND. verbose) THEN                   ! Control print
1394         WRITE(numout,*) 
1395         WRITE(numout,*) 'dia_wri_region_mean finished for ', trim(tmp_name)
1396         WRITE(numout,*) 
1397         CALL FLUSH(numout)
1398      ENDIF
1399     
1400   END SUBROUTINE dia_wri_region_mean
1401
1402
1403   !!======================================================================
1404END MODULE diaregmean
Note: See TracBrowser for help on using the repository browser.