New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diaregmean.F90 in branches/UKMO/CO6_shelfclimate_fabm/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/CO6_shelfclimate_fabm/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90 @ 10493

Last change on this file since 10493 was 10493, checked in by hadjt, 5 years ago

Updated region mean to fix processor overlap double count error

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