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

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

Region means
IOM_use to only process variables in use.
namelist parameter to control whether hourly means, or instanteous values
Instantaneous values of heat, salt mass and volume.

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