source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90 @ 11038

Last change on this file since 11038 was 11038, checked in by rrenshaw, 17 months ago

code fixes

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