source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis3/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90 @ 11070

Last change on this file since 11070 was 11070, checked in by rrenshaw, 18 months ago

code fixes

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