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

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

DIA/diaregmean.F90

Region mean namelist switch added to allow area weighted region means.

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