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

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

Region mean.
fixing the total and count dianostic output, and adding a total x area region mean diagnostic

File size: 66.8 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_totarea,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,area_mat,totarea_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_totarea(n_regions_output),  STAT= ierr )
945        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_totarea array' )
946      ALLOCATE( zrmet_var(n_regions_output),  STAT= ierr )
947        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_var array' )
948      ALLOCATE( zrmet_cnt(n_regions_output),  STAT= ierr )
949        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_cnt array' )
950      ALLOCATE( zrmet_mask_id(n_regions_output),  STAT= ierr )
951        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_mask_id array' )
952      ALLOCATE( zrmet_reg_id(n_regions_output),  STAT= ierr )
953        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
954
955
956      ALLOCATE( zrmet_min(n_regions_output),  STAT= ierr )
957        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_min array' )
958      ALLOCATE( zrmet_max(n_regions_output),  STAT= ierr )
959        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_max array' )
960     
961      ALLOCATE( zrmet_out(n_regions_output),  STAT= ierr )
962        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
963
964 
965     
966        IF(lwp .AND. verbose) THEN
967              WRITE(numout,*)
968              WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//';'
969              WRITE(numout,*)
970        ENDIF
971       
972      DO ji = 1,jpi
973          DO jj = 1,jpj
974            internal_infield(ji,jj) = infield(ji,jj)
975          END DO
976      END DO
977       
978      ! Check for NANS # JT   03/09/2018
979      DO ji = 1,jpi
980          DO jj = 1,jpj
981              IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
982                  IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
983                      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)
984                      internal_infield(ji,jj) = 0.
985                  ENDIF
986              ELSE               
987                  IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
988                      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)
989                      internal_infield(ji,jj) = 0.
990                  ENDIF
991              ENDIF
992          END DO
993      END DO
994     
995     
996      zrmet_ave(:) = zmdi
997      zrmet_tot(:) = zmdi
998      zrmet_totarea(:) = zmdi
999      zrmet_var(:) = zmdi
1000      zrmet_cnt(:) = zmdi
1001      zrmet_mask_id(:) = zmdi
1002      zrmet_reg_id(:) = zmdi
1003     
1004      zrmet_min(:) = zmdi
1005      zrmet_max(:) = zmdi
1006      reg_ind_cnt = 1
1007     
1008     
1009      ! loop though the masks
1010      DO maskno = 1,nmasks
1011          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin mask loops: ',maskno
1012         
1013         
1014          ! For each mask, get the number of regions (nreg), and a local copy of the region.
1015          nreg = nreg_mat(maskno)
1016          internal_region_mask = region_mask(:,:,maskno)
1017         
1018          ! allocate temporary stat arrays, and set to zero
1019          ALLOCATE( ave_mat(nreg),  STAT= ierr )
1020          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ave_mat array' )
1021          ALLOCATE( tot_mat(nreg),      STAT= ierr )
1022          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate tot_mat array' )
1023          ALLOCATE( num_mat(nreg),  STAT= ierr )
1024          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate num_mat array' )
1025          ALLOCATE( var_mat(nreg),  STAT= ierr )
1026          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate var_mat array' )
1027          ALLOCATE( ssq_mat(nreg),  STAT= ierr )
1028          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ssq_mat array' )
1029          ALLOCATE( cnt_mat(nreg),  STAT= ierr )
1030          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate cnt_mat array' )
1031          ALLOCATE( area_mat(nreg),  STAT= ierr )
1032          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate area_mat array' )
1033          ALLOCATE( totarea_mat(nreg),  STAT= ierr )
1034          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate totarea_mat array' )
1035
1036          !ALLOCATE( min_mat(nreg),  STAT= ierr )
1037          !IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate min_mat array' )
1038          !ALLOCATE( max_mat(nreg),  STAT= ierr )
1039          !IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate max_mat array' )
1040
1041          ALLOCATE( reg_id_mat(nreg),  STAT= ierr )
1042          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate reg_id_mat array' )
1043          ALLOCATE( mask_id_mat(nreg),  STAT= ierr )
1044          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate mask_id_mat array' )
1045
1046
1047         
1048          ave_mat(:) = 0.
1049          tot_mat(:) = 0.
1050          num_mat(:) = 0.
1051          var_mat(:) = 0.
1052          cnt_mat(:) = 0.
1053          ssq_mat(:) = 0.
1054          area_mat(:) = 0.
1055          totarea_mat(:) = 0.
1056
1057
1058          !min_mat(:) = zmdi
1059          !max_mat(:) = -zmdi
1060          reg_id_mat(:) = 0.
1061          mask_id_mat(:) = 0.
1062         
1063          ! loop though the array. for each sea grid box where tmask == 1),
1064          ! read which region the grid box is in, add the value of the gridbox (and its square)
1065          ! to the total for that region, and then increment the counter for that region.
1066          !CALL cpu_time(start_reg_mean_loop)
1067          !WRITE(numout,*) kt,start_reg_mean_loop
1068          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin spatial loops: '
1069          DO ji = nldi,nlei
1070              DO jj = nldj,nlej
1071                    IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
1072                        ind = internal_region_mask(ji,jj)+1
1073                        !tot_mat(ind) = tot_mat(ind) + (internal_infield(ji,jj))
1074                        !ssq_mat(ind) = ssq_mat(ind) + ( internal_infield(ji,jj) *  internal_infield(ji,jj))
1075                        !cnt_mat(ind) = cnt_mat(ind) + 1.
1076                        ! Area Weighted values - region_area_mat == 1. or area depending on ln_diaregmean_areawgt
1077                        totarea_mat(ind) = totarea_mat(ind) + (region_area_mat(ji,jj)*internal_infield(ji,jj))
1078                        tot_mat(ind) = tot_mat(ind) + (internal_infield(ji,jj))
1079                        ssq_mat(ind) = ssq_mat(ind) + (region_area_mat(ji,jj)*(internal_infield(ji,jj) * internal_infield(ji,jj)))
1080                        cnt_mat(ind) = cnt_mat(ind) + 1.
1081                        area_mat(ind) = area_mat(ind) + (region_area_mat(ji,jj)*1.)
1082
1083
1084
1085                        !min_mat(ind) = min(min_mat(ind),internal_infield(ji,jj))
1086                        !max_mat(ind) = max(max_mat(ind),internal_infield(ji,jj))
1087                    ENDIF
1088              END DO
1089          END DO
1090          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finish spatial loops: '
1091          ! sum the totals, the counts, and the squares across the processors         
1092          CALL mpp_sum( 'diaregionmean',tot_mat,nreg )
1093          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum tot'
1094          CALL mpp_sum( 'diaregionmean',cnt_mat,nreg )
1095          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum cnt'
1096          CALL mpp_sum( 'diaregionmean',area_mat,nreg )
1097          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum area'
1098          CALL mpp_sum( 'diaregionmean',totarea_mat,nreg )
1099          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum totarea_mat'
1100
1101
1102
1103          !tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_var'))
1104          !IF (iom_use(trim(tmp_name_iom)) .OR. ln_diaregmean_bin .OR. ln_diaregmean_ascii) THEN
1105              CALL mpp_sum( 'diaregionmean',ssq_mat,nreg )
1106              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum ssq'
1107          !ENDIF
1108         
1109   
1110
1111          !tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_min'))
1112          !IF (iom_use(trim(tmp_name_iom)) .OR. ln_diaregmean_bin .OR. ln_diaregmean_ascii) THEN
1113              !CALL mpp_min( 'diaregionmean',min_mat,nreg )
1114              !IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_min'
1115          !ENDIF
1116         
1117
1118          !tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_max'))
1119          !IF (iom_use(trim(tmp_name_iom)) .OR. ln_diaregmean_bin .OR. ln_diaregmean_ascii)  THEN
1120              !CALL mpp_max( 'diaregionmean',max_mat,nreg )
1121              !IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_max'
1122          !ENDIF
1123         
1124          !calculate the mean and variance from the total, sum of squares and the count.
1125          ! When area weighting, you can't area weight the total.
1126          ! this if block may be redundant, as totarea_mat == tot_mat, and cnt_mat == area_mat when ln_diaregmean_areawgt == False
1127          IF (ln_diaregmean_areawgt) THEN
1128            ave_mat = totarea_mat(:)/area_mat(:)
1129            var_mat = ssq_mat(:)/area_mat(:) - (ave_mat(:)*ave_mat(:))
1130          ELSE
1131            ave_mat = tot_mat(:)/cnt_mat(:)
1132            var_mat = ssq_mat(:)/cnt_mat(:) - (ave_mat(:)*ave_mat(:))
1133          ENDIF
1134         
1135         
1136         
1137          !mask array of mask and region number.
1138          DO jj = 1,nreg
1139              reg_id_mat(jj) = real(jj-1)
1140              mask_id_mat(jj) = real(maskno)
1141          END DO
1142         
1143         
1144          !write text and binary, and note region statistics for current mask for later iom_put
1145          IF( lwp ) THEN 
1146         
1147              !Write out ascii and binary if requred
1148              IF ( ln_diaregmean_bin ) THEN
1149                  !Writing out regional mean time series to binary files
1150                  WRITE(numdct_reg_bin) tmp_name,kt,maskno,n_regions_output
1151                  WRITE(numdct_reg_bin) ave_mat
1152                  WRITE(numdct_reg_bin) tot_mat
1153                  WRITE(numdct_reg_bin) var_mat
1154                  WRITE(numdct_reg_bin) ssq_mat
1155                  WRITE(numdct_reg_bin) cnt_mat
1156                  !WRITE(numdct_reg_bin) min_mat
1157                  !WRITE(numdct_reg_bin) max_mat
1158              ENDIF
1159             
1160              IF ( ln_diaregmean_ascii  ) THEN
1161                  !Writing out regional mean time series to text files
1162
1163                  WRITE(nreg_string, "(I5)") nreg
1164                  FormatString = "(A30,"//trim(nreg_string)//"F25.3)"
1165                  WRITE(numdct_reg_txt, FMT="(A30,I6,I6)") tmp_name,kt,maskno           
1166                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ave_mat:", ave_mat
1167                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"tot_mat:", tot_mat
1168                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"var_mat:", var_mat
1169                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ssq_mat:", ssq_mat
1170                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"cnt_mat:", cnt_mat
1171                  !WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"min_mat:", min_mat
1172                  !WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"max_mat:", max_mat
1173                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"reg_mat:", reg_id_mat
1174                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"msk_mat:", mask_id_mat
1175
1176              ENDIF
1177          ENDIF
1178             
1179          ! JT Fixed, was not meant to be inside the lwp if block
1180          DO jm = 1,nreg
1181              zrmet_ave(    reg_ind_cnt) =     ave_mat(jm)
1182              zrmet_tot(    reg_ind_cnt) =     tot_mat(jm)
1183              zrmet_totarea(    reg_ind_cnt) =     totarea_mat(jm)
1184              zrmet_var(    reg_ind_cnt) =     var_mat(jm)
1185              zrmet_cnt(    reg_ind_cnt) =     cnt_mat(jm)
1186              !zrmet_min(    reg_ind_cnt) =     min_mat(jm)
1187              !zrmet_max(    reg_ind_cnt) =     max_mat(jm)
1188              zrmet_reg_id( reg_ind_cnt) =  reg_id_mat(jm)
1189              zrmet_mask_id(reg_ind_cnt) = mask_id_mat(jm)
1190           
1191              reg_ind_cnt = reg_ind_cnt + 1 
1192          END DO
1193     
1194       
1195          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean about to deallocated arrays for ',kt,maskno
1196          DEALLOCATE(ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,reg_id_mat,mask_id_mat,totarea_mat, area_mat)
1197          !DEALLOCATE(min_mat,max_mat)
1198
1199          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean deallocated arrays for ',kt,maskno
1200          IF(lwp .AND. ln_diaregmean_ascii ) CALL FLUSH(numdct_reg_txt)
1201          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean flushed region mean text for ',kt,maskno
1202      END DO
1203
1204      IF(lwp .AND. verbose) THEN                   ! Control print
1205         WRITE(numout,*) 'dia_regmean ready to start iom_put'
1206         CALL FLUSH(numout)
1207      ENDIF
1208     
1209      !With current field_def.xml and iodef.xml, these fields must be output, so set to dummy values if not required.
1210     
1211      IF ( ln_diaregmean_nc ) THEN
1212     
1213          zrmet_out(:) = 0
1214          zrmet_val = 0
1215          tmp_name_iom = ''
1216
1217          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean ready to start iom_put: ',trim(tmp_name)
1218         
1219          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_ave'))
1220          IF (iom_use(trim(tmp_name_iom))) THEN
1221              DO jm = 1,n_regions_output
1222                zrmet_val = zrmet_ave(jm)
1223    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1224    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1225                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1226                zrmet_out(jm) = zrmet_val
1227              END DO     
1228              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1229              CALL iom_put(trim(tmp_name_iom), zrmet_out(:) ) 
1230              zrmet_out(:) = 0
1231              zrmet_val = 0
1232              tmp_name_iom = ''
1233          ENDIF
1234
1235          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_tot'))
1236          IF (iom_use(trim(tmp_name_iom))) THEN
1237              DO jm = 1,n_regions_output
1238                zrmet_val = zrmet_tot(jm)
1239    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1240    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1241                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1242                zrmet_out(jm) = zrmet_val
1243              END DO
1244              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1245              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) ) 
1246              zrmet_out(:) = 0
1247              zrmet_val = 0
1248              tmp_name_iom = ''
1249          ENDIF
1250
1251
1252
1253          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_totarea'))
1254          IF (iom_use(trim(tmp_name_iom))) THEN
1255              DO jm = 1,n_regions_output
1256                zrmet_val = zrmet_totarea(jm)
1257    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1258    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1259                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1260                zrmet_out(jm) = zrmet_val
1261              END DO
1262              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1263              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) ) 
1264              zrmet_out(:) = 0
1265              zrmet_val = 0
1266              tmp_name_iom = ''
1267          ENDIF
1268
1269          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_var'))
1270          IF (iom_use(trim(tmp_name_iom))) THEN
1271              DO jm = 1,n_regions_output
1272                zrmet_val = zrmet_var(jm)
1273    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1274    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1275                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1276                zrmet_out(jm) = zrmet_val
1277              END DO
1278              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1279              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) )
1280              zrmet_out(:) = 0
1281              zrmet_val = 0
1282              tmp_name_iom = ''
1283          ENDIF
1284
1285          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))
1286          IF (iom_use(trim(tmp_name_iom))) THEN
1287              DO jm = 1,n_regions_output
1288                zrmet_val = zrmet_cnt(jm)
1289    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1290    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1291                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1292                zrmet_out(jm) = zrmet_val
1293              END DO
1294              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1295              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) )
1296              zrmet_out(:) = 0
1297              zrmet_val = 0
1298              tmp_name_iom = ''
1299          ENDIF
1300
1301          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_reg_id'))
1302          IF (iom_use(trim(tmp_name_iom))) THEN
1303              DO jm = 1,n_regions_output
1304                zrmet_val = zrmet_reg_id(jm)
1305    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1306    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1307                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1308                zrmet_out(jm) = zrmet_val
1309              END DO
1310              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom), zrmet_out(1)
1311              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) ) 
1312              zrmet_out(:) = 0
1313              zrmet_val = 0
1314              tmp_name_iom = ''
1315          ENDIF
1316
1317          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_mask_id'))
1318          IF (iom_use(trim(tmp_name_iom))) THEN
1319              DO jm = 1,n_regions_output
1320                zrmet_val = zrmet_mask_id(jm)
1321    !            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
1322    !            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
1323                if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
1324                zrmet_out(jm) = zrmet_val
1325              END DO
1326              IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)         
1327              CALL iom_put( trim(tmp_name_iom), zrmet_out(:) )
1328              zrmet_out(:) = 0
1329              zrmet_val = 0
1330              tmp_name_iom = ''
1331          ENDIF
1332     
1333!      ELSE
1334!       
1335!          ALLOCATE( dummy_zrmet(jpi,jpj,n_regions_output),  STAT= ierr )
1336!            IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate dummy_zrmet array' )
1337
1338!          DO jm = 1,n_regions_output
1339!              dummy_zrmet(:,:,jm) =     real(jm,wp)
1340!          END DO
1341
1342!          DO jm = 1,9
1343!              !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 )
1344!              !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 )
1345!              !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 )
1346!              !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 )
1347!              !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 )
1348!              !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 )
1349
1350
1351!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_ave'))
1352!              IF (iom_use(trim(tmp_name_iom))) THEN
1353!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1354!              ENDIF
1355!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_tot'))
1356!              IF (iom_use(trim(tmp_name_iom))) THEN
1357!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1358!              ENDIF
1359!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_var'))
1360!              IF (iom_use(trim(tmp_name_iom))) THEN
1361!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1362!              ENDIF
1363!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))
1364!              IF (iom_use(trim(tmp_name_iom))) THEN
1365!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1366!              ENDIF
1367!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_reg_id'))
1368!              IF (iom_use(trim(tmp_name_iom))) THEN
1369!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1370!              ENDIF
1371!              tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_mask_id'))
1372!              IF (iom_use(trim(tmp_name_iom))) THEN
1373!                 CALL iom_put( trim(tmp_name_iom), dummy_zrmet(1,1,:) ) !dummy_zrmet(1,1,:) ) )
1374!              ENDIF
1375!
1376!          END DO
1377!   
1378!          DEALLOCATE( dummy_zrmet)
1379      ENDIF
1380     
1381      DEALLOCATE(zrmet_ave,zrmet_tot,zrmet_totarea,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id,zrmet_min,zrmet_max,zrmet_out)
1382
1383      IF(lwp .AND. verbose) THEN                   ! Control print
1384         WRITE(numout,*) 
1385         WRITE(numout,*) 'dia_wri_region_mean finished for ', trim(tmp_name)
1386         WRITE(numout,*) 
1387         CALL FLUSH(numout)
1388      ENDIF
1389     
1390   END SUBROUTINE dia_wri_region_mean
1391
1392
1393   !!======================================================================
1394END MODULE diaregmean
Note: See TracBrowser for help on using the repository browser.