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.
diaptr.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 30.9 KB
Line 
1MODULE diaptr
2   !!======================================================================
3   !!                       ***  MODULE  diaptr  ***
4   !! Ocean physics:  Computes meridonal transports and zonal means
5   !!=====================================================================
6   !! History :  1.0  ! 2003-09  (C. Talandier, G. Madec)  Original code
7   !!            2.0  ! 2006-01  (A. Biastoch)  Allow sub-basins computation
8   !!            3.2  ! 2010-03  (O. Marti, S. Flavoni) Add fields
9   !!            3.3  ! 2010-10  (G. Madec)  dynamical allocation
10   !!            3.6  ! 2014-12  (C. Ethe) use of IOM
11   !!            3.6  ! 2016-06  (T. Graham) Addition of diagnostics for CMIP6
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   dia_ptr      : Poleward Transport Diagnostics module
16   !!   dia_ptr_init : Initialization, namelist read
17   !!   ptr_sjk      : "zonal" mean computation of a field - tracer or flux array
18   !!   ptr_sj       : "zonal" and vertical sum computation of a "meridional" flux array
19   !!                   (Generic interface to ptr_sj_3d, ptr_sj_2d)
20   !!----------------------------------------------------------------------
21   USE oce              ! ocean dynamics and active tracers
22   USE dom_oce          ! ocean space and time domain
23   USE phycst           ! physical constants
24   !
25   USE iom              ! IOM library
26   USE in_out_manager   ! I/O manager
27   USE lib_mpp          ! MPP library
28   USE timing           ! preformance summary
29
30   IMPLICIT NONE
31   PRIVATE
32
33   INTERFACE ptr_sj
34      MODULE PROCEDURE ptr_sj_3d, ptr_sj_2d
35   END INTERFACE
36
37   PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines
38   PUBLIC   ptr_sjk        !
39   PUBLIC   dia_ptr_init   ! call in step module
40   PUBLIC   dia_ptr        ! call in step module
41   PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines
42
43   !                                  !!** namelist  namptr  **
44   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_adv, htr_ldf, htr_eiv   !: Heat TRansports (adv, diff, Bolus.)
45   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   str_adv, str_ldf, str_eiv   !: Salt TRansports (adv, diff, Bolus.)
46   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_ove, str_ove   !: heat Salt TRansports ( overturn.)
47   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_btr, str_btr   !: heat Salt TRansports ( barotropic )
48
49   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F)
50   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation
51   INTEGER, PUBLIC ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T)
52
53   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
54   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rau0 x Cp)
55   REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg
56
57   CHARACTER(len=3), ALLOCATABLE, SAVE, DIMENSION(:)     :: clsubb
58   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
59   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   :: btm30   ! mask out Southern Ocean (=0 south of 30°S)
60
61   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)     :: p_fval1d
62   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: p_fval2d
63
64   !! * Substitutions
65#  include "vectopt_loop_substitute.h90"
66   !!----------------------------------------------------------------------
67   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
68   !! $Id$
69   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
70   !!----------------------------------------------------------------------
71CONTAINS
72
73   SUBROUTINE dia_ptr( pvtr )
74      !!----------------------------------------------------------------------
75      !!                  ***  ROUTINE dia_ptr  ***
76      !!----------------------------------------------------------------------
77      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
78      !
79      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
80      REAL(wp) ::   zsfc,zvfc               ! local scalar
81      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
82      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace
83      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace
84      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace
85      REAL(wp), DIMENSION(jpj)     ::  vsum   ! 1D workspace
86      REAL(wp), DIMENSION(jpj,jpts)     ::  tssum   ! 1D workspace
87 
88      !
89      !overturning calculation
90      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   sjk  , r1_sjk ! i-mean i-k-surface and its inverse
91      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   v_msf, sn_jk  , tn_jk ! i-mean T and S, j-Stream-Function
92      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace
93
94
95      CHARACTER( len = 12 )  :: cl1
96      !!----------------------------------------------------------------------
97      !
98      IF( nn_timing == 1 )   CALL timing_start('dia_ptr')
99
100      !
101      IF( PRESENT( pvtr ) ) THEN
102         IF( iom_use("zomsfglo") ) THEN    ! effective MSF
103            z3d(1,:,:) = ptr_sjk( pvtr(:,:,:) )  ! zonal cumulative effective transport
104            DO jk = 2, jpkm1 
105              z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)   ! effective j-Stream-Function (MSF)
106            END DO
107            DO ji = 1, jpi
108               z3d(ji,:,:) = z3d(1,:,:)
109            ENDDO
110            cl1 = TRIM('zomsf'//clsubb(1) )
111            CALL iom_put( cl1, z3d * rc_sv )
112            DO jn = 2, nptr                                    ! by sub-basins
113               z3d(1,:,:) =  ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) 
114               DO jk = 2, jpkm1 
115                  z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)    ! effective j-Stream-Function (MSF)
116               END DO
117               DO ji = 1, jpi
118                  z3d(ji,:,:) = z3d(1,:,:)
119               ENDDO
120               cl1 = TRIM('zomsf'//clsubb(jn) )
121               CALL iom_put( cl1, z3d * rc_sv )
122            END DO
123         ENDIF
124         IF( iom_use("sopstove") .OR. iom_use("sophtove") .OR. iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN
125            ! define fields multiplied by scalar
126            zmask(:,:,:) = 0._wp
127            zts(:,:,:,:) = 0._wp
128            zvn(:,:,:) = 0._wp
129            DO jk = 1, jpkm1
130               DO jj = 1, jpjm1
131                  DO ji = 1, jpi
132                     zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk)
133                     zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
134                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
135                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
136                     zvn(ji,jj,jk)        = vn(ji,jj,jk)         * zvfc
137                  ENDDO
138               ENDDO
139             ENDDO
140         ENDIF
141         IF( iom_use("sopstove") .OR. iom_use("sophtove") ) THEN
142             sjk(:,:,1) = ptr_sjk( zmask(:,:,:), btmsk(:,:,1) )
143             r1_sjk(:,:,1) = 0._wp
144             WHERE( sjk(:,:,1) /= 0._wp )   r1_sjk(:,:,1) = 1._wp / sjk(:,:,1)
145
146             ! i-mean T and S, j-Stream-Function, global
147             tn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_tem) ) * r1_sjk(:,:,1)
148             sn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_sal) ) * r1_sjk(:,:,1)
149             v_msf(:,:,1) = ptr_sjk( zvn(:,:,:) )
150
151             htr_ove(:,1) = SUM( v_msf(:,:,1)*tn_jk(:,:,1) ,2 )
152             str_ove(:,1) = SUM( v_msf(:,:,1)*sn_jk(:,:,1) ,2 )
153
154             z2d(1,:) = htr_ove(:,1) * rc_pwatt        !  (conversion in PW)
155             DO ji = 1, jpi
156               z2d(ji,:) = z2d(1,:)
157             ENDDO
158             cl1 = 'sophtove'
159             CALL iom_put( TRIM(cl1), z2d )
160             z2d(1,:) = str_ove(:,1) * rc_ggram        !  (conversion in Gg)
161             DO ji = 1, jpi
162               z2d(ji,:) = z2d(1,:)
163             ENDDO
164             cl1 = 'sopstove'
165             CALL iom_put( TRIM(cl1), z2d )
166             IF( ln_subbas ) THEN
167                DO jn = 2, nptr
168                    sjk(:,:,jn) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
169                    r1_sjk(:,:,jn) = 0._wp
170                    WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
171
172                    ! i-mean T and S, j-Stream-Function, basin
173                    tn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
174                    sn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
175                    v_msf(:,:,jn) = ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) ) 
176                    htr_ove(:,jn) = SUM( v_msf(:,:,jn)*tn_jk(:,:,jn) ,2 )
177                    str_ove(:,jn) = SUM( v_msf(:,:,jn)*sn_jk(:,:,jn) ,2 )
178
179                    z2d(1,:) = htr_ove(:,jn) * rc_pwatt !  (conversion in PW)
180                    DO ji = 1, jpi
181                        z2d(ji,:) = z2d(1,:)
182                    ENDDO
183                    cl1 = TRIM('sophtove_'//clsubb(jn))
184                    CALL iom_put( cl1, z2d )
185                    z2d(1,:) = str_ove(:,jn) * rc_ggram        ! (conversion in Gg)
186                    DO ji = 1, jpi
187                        z2d(ji,:) = z2d(1,:)
188                    ENDDO
189                    cl1 = TRIM('sopstove_'//clsubb(jn))
190                    CALL iom_put( cl1, z2d )
191                END DO
192             ENDIF
193         ENDIF
194         IF( iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN
195         ! Calculate barotropic heat and salt transport here
196             sjk(:,1,1) = ptr_sj( zmask(:,:,:), btmsk(:,:,1) )
197             r1_sjk(:,1,1) = 0._wp
198             WHERE( sjk(:,1,1) /= 0._wp )   r1_sjk(:,1,1) = 1._wp / sjk(:,1,1)
199           
200            vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,1))
201            tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,1) )
202            tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,1) )
203            htr_btr(:,1) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,1)
204            str_btr(:,1) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,1)
205            z2d(1,:) = htr_btr(:,1) * rc_pwatt        !  (conversion in PW)
206            DO ji = 2, jpi
207               z2d(ji,:) = z2d(1,:)
208            ENDDO
209            cl1 = 'sophtbtr'
210            CALL iom_put( TRIM(cl1), z2d )
211            z2d(1,:) = str_btr(:,1) * rc_ggram        !  (conversion in Gg)
212            DO ji = 2, jpi
213              z2d(ji,:) = z2d(1,:)
214            ENDDO
215            cl1 = 'sopstbtr'
216            CALL iom_put( TRIM(cl1), z2d )
217            IF( ln_subbas ) THEN
218                DO jn = 2, nptr
219                    sjk(:,1,jn) = ptr_sj( zmask(:,:,:), btmsk(:,:,jn) )
220                    r1_sjk(:,1,jn) = 0._wp
221                    WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
222                    vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,jn))
223                    tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
224                    tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
225                    htr_btr(:,jn) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,jn)
226                    str_btr(:,jn) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,jn)
227                    z2d(1,:) = htr_btr(:,jn) * rc_pwatt !  (conversion in PW)
228                    DO ji = 1, jpi
229                        z2d(ji,:) = z2d(1,:)
230                    ENDDO
231                    cl1 = TRIM('sophtbtr_'//clsubb(jn))
232                    CALL iom_put( cl1, z2d )
233                    z2d(1,:) = str_btr(:,jn) * rc_ggram        ! (conversion in Gg)
234                    DO ji = 1, jpi
235                        z2d(ji,:) = z2d(1,:)
236                    ENDDO
237                    cl1 = TRIM('sopstbtr_'//clsubb(jn))
238                    CALL iom_put( cl1, z2d )
239               ENDDO
240            ENDIF !ln_subbas
241         ENDIF !iom_use("sopstbtr....)
242         !
243      ELSE
244         !
245         IF( iom_use("zotemglo") ) THEN    ! i-mean i-k-surface
246            DO jk = 1, jpkm1
247               DO jj = 1, jpj
248                  DO ji = 1, jpi
249                     zsfc = e1t(ji,jj) * e3t_n(ji,jj,jk)
250                     zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
251                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc
252                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc
253                  END DO
254               END DO
255            END DO
256            DO jn = 1, nptr
257               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
258               cl1 = TRIM('zosrf'//clsubb(jn) )
259               CALL iom_put( cl1, zmask )
260               !
261               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) &
262                  &            / MAX( zmask(1,:,:), 10.e-15 )
263               DO ji = 1, jpi
264                  z3d(ji,:,:) = z3d(1,:,:)
265               ENDDO
266               cl1 = TRIM('zotem'//clsubb(jn) )
267               CALL iom_put( cl1, z3d )
268               !
269               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) &
270                  &            / MAX( zmask(1,:,:), 10.e-15 )
271               DO ji = 1, jpi
272                  z3d(ji,:,:) = z3d(1,:,:)
273               ENDDO
274               cl1 = TRIM('zosal'//clsubb(jn) )
275               CALL iom_put( cl1, z3d )
276            END DO
277         ENDIF
278         !
279         !                                ! Advective and diffusive heat and salt transport
280         IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN   
281            z2d(1,:) = htr_adv(:,1) * rc_pwatt        !  (conversion in PW)
282            DO ji = 1, jpi
283               z2d(ji,:) = z2d(1,:)
284            ENDDO
285            cl1 = 'sophtadv'                 
286            CALL iom_put( TRIM(cl1), z2d )
287            z2d(1,:) = str_adv(:,1) * rc_ggram        ! (conversion in Gg)
288            DO ji = 1, jpi
289               z2d(ji,:) = z2d(1,:)
290            ENDDO
291            cl1 = 'sopstadv'
292            CALL iom_put( TRIM(cl1), z2d )
293            IF( ln_subbas ) THEN
294              DO jn=2,nptr
295               z2d(1,:) = htr_adv(:,jn) * rc_pwatt        !  (conversion in PW)
296               DO ji = 1, jpi
297                 z2d(ji,:) = z2d(1,:)
298               ENDDO
299               cl1 = TRIM('sophtadv_'//clsubb(jn))                 
300               CALL iom_put( cl1, z2d )
301               z2d(1,:) = str_adv(:,jn) * rc_ggram        ! (conversion in Gg)
302               DO ji = 1, jpi
303                  z2d(ji,:) = z2d(1,:)
304               ENDDO
305               cl1 = TRIM('sopstadv_'//clsubb(jn))                 
306               CALL iom_put( cl1, z2d )             
307              ENDDO
308            ENDIF
309         ENDIF
310         !
311         IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN   
312            z2d(1,:) = htr_ldf(:,1) * rc_pwatt        !  (conversion in PW)
313            DO ji = 1, jpi
314               z2d(ji,:) = z2d(1,:)
315            ENDDO
316            cl1 = 'sophtldf'
317            CALL iom_put( TRIM(cl1), z2d )
318            z2d(1,:) = str_ldf(:,1) * rc_ggram        !  (conversion in Gg)
319            DO ji = 1, jpi
320               z2d(ji,:) = z2d(1,:)
321            ENDDO
322            cl1 = 'sopstldf'
323            CALL iom_put( TRIM(cl1), z2d )
324            IF( ln_subbas ) THEN
325              DO jn=2,nptr
326               z2d(1,:) = htr_ldf(:,jn) * rc_pwatt        !  (conversion in PW)
327               DO ji = 1, jpi
328                 z2d(ji,:) = z2d(1,:)
329               ENDDO
330               cl1 = TRIM('sophtldf_'//clsubb(jn))                 
331               CALL iom_put( cl1, z2d )
332               z2d(1,:) = str_ldf(:,jn) * rc_ggram        ! (conversion in Gg)
333               DO ji = 1, jpi
334                  z2d(ji,:) = z2d(1,:)
335               ENDDO
336               cl1 = TRIM('sopstldf_'//clsubb(jn))                 
337               CALL iom_put( cl1, z2d )             
338              ENDDO
339            ENDIF
340         ENDIF
341
342         IF( iom_use("sophteiv") .OR. iom_use("sopsteiv") ) THEN
343            z2d(1,:) = htr_eiv(:,1) * rc_pwatt        !  (conversion in PW)
344            DO ji = 1, jpi
345               z2d(ji,:) = z2d(1,:)
346            ENDDO
347            cl1 = 'sophteiv'
348            CALL iom_put( TRIM(cl1), z2d )
349            z2d(1,:) = str_eiv(:,1) * rc_ggram        !  (conversion in Gg)
350            DO ji = 1, jpi
351               z2d(ji,:) = z2d(1,:)
352            ENDDO
353            cl1 = 'sopsteiv'
354            CALL iom_put( TRIM(cl1), z2d )
355            IF( ln_subbas ) THEN
356               DO jn=2,nptr
357                  z2d(1,:) = htr_eiv(:,jn) * rc_pwatt        !  (conversion in PW)
358                  DO ji = 1, jpi
359                     z2d(ji,:) = z2d(1,:)
360                  ENDDO
361                  cl1 = TRIM('sophteiv_'//clsubb(jn))                 
362                  CALL iom_put( cl1, z2d )
363                  z2d(1,:) = str_eiv(:,jn) * rc_ggram        ! (conversion in Gg)
364                  DO ji = 1, jpi
365                     z2d(ji,:) = z2d(1,:)
366                  ENDDO
367                  cl1 = TRIM('sopsteiv_'//clsubb(jn)) 
368                  CALL iom_put( cl1, z2d )             
369               ENDDO
370            ENDIF
371         ENDIF
372         !
373      ENDIF
374      !
375      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr')
376      !
377   END SUBROUTINE dia_ptr
378
379
380   SUBROUTINE dia_ptr_init
381      !!----------------------------------------------------------------------
382      !!                  ***  ROUTINE dia_ptr_init  ***
383      !!                   
384      !! ** Purpose :   Initialization, namelist read
385      !!----------------------------------------------------------------------
386      INTEGER ::  jn, jj, ji   ! local integers
387      INTEGER ::  inum, ierr   ! local integers
388      INTEGER ::  ios          ! Local integer output status for namelist read
389      !!
390      NAMELIST/namptr/ ln_diaptr, ln_subbas
391      !!----------------------------------------------------------------------
392
393      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
394      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901)
395901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist', lwp )
396
397      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport
398      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 )
399902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp )
400      IF(lwm) WRITE ( numond, namptr )
401
402      IF(lwp) THEN                     ! Control print
403         WRITE(numout,*)
404         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
405         WRITE(numout,*) '~~~~~~~~~~~~'
406         WRITE(numout,*) '   Namelist namptr : set ptr parameters'
407         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr
408         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas
409      ENDIF
410
411      IF( ln_diaptr ) THEN 
412         !
413         IF( ln_subbas ) THEN
414            nptr = 5            ! Global, Atlantic, Pacific, Indian, Indo-Pacific
415            ALLOCATE( clsubb(nptr) )
416            clsubb(1) = 'glo' ;  clsubb(2) = 'atl'  ;  clsubb(3) = 'pac'  ;  clsubb(4) = 'ind'  ;  clsubb(5) = 'ipc'
417         ELSE               
418            nptr = 1       ! Global only
419            ALLOCATE( clsubb(nptr) )
420            clsubb(1) = 'glo' 
421         ENDIF
422
423         !                                      ! allocate dia_ptr arrays
424         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
425
426         rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt
427
428         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
429
430         IF( ln_subbas ) THEN                ! load sub-basin mask
431            CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  )
432            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
433            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
434            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
435            CALL iom_close( inum )
436!$OMP PARALLEL DO schedule(static) private(jj,ji)
437            DO jj = 1, jpj
438               DO ji = 1, jpi
439                  btmsk(ji,jj,5) = MAX ( btmsk(ji,jj,3), btmsk(ji,jj,4) )          ! Indo-Pacific basin
440                  IF( gphit(ji,jj) < -30._wp) THEN   ;   btm30(ji,jj) = 0._wp      ! mask out Southern Ocean
441                  ELSE                               ;   btm30(ji,jj) = ssmask(ji,jj)
442                  END IF
443               END DO
444            END DO
445         ENDIF
446   
447!$OMP PARALLEL
448!$OMP DO schedule(static) private(jj,ji)
449         DO jj = 1, jpj
450           DO ji = 1, jpi
451              btmsk(ji,jj,1) = tmask_i(ji,jj)                          ! global ocean
452           END DO
453         END DO
454     
455         DO jn = 1, nptr
456!$OMP DO schedule(static) private(jj,ji)
457            DO jj = 1, jpj
458               DO ji = 1, jpi
459                  btmsk(ji,jj,jn) = btmsk(ji,jj,jn) * tmask_i(ji,jj)               ! interior domain only
460               END DO
461            END DO
462         END DO
463
464         ! Initialise arrays to zero because diatpr is called before they are first calculated
465         ! Note that this means diagnostics will not be exactly correct when model run is restarted.
466!$OMP DO schedule(static) private(jj,ji)
467         DO jj = 1, jpj
468            DO ji = 1, jpi
469               htr_adv(ji,jj) = 0._wp  ;  str_adv(ji,jj) =  0._wp 
470               htr_ldf(ji,jj) = 0._wp  ;  str_ldf(ji,jj) =  0._wp 
471               htr_eiv(ji,jj) = 0._wp  ;  str_eiv(ji,jj) =  0._wp 
472               htr_ove(ji,jj) = 0._wp  ;   str_ove(ji,jj) =  0._wp
473               htr_btr(ji,jj) = 0._wp  ;   str_btr(ji,jj) =  0._wp
474             END DO
475         END DO
476              !
477!$OMP END PARALLEL
478      ENDIF 
479      !
480   END SUBROUTINE dia_ptr_init
481
482   SUBROUTINE dia_ptr_hst( ktra, cptr, pva ) 
483      !!----------------------------------------------------------------------
484      !!                    ***  ROUTINE dia_ptr_hst ***
485      !!----------------------------------------------------------------------
486      !! Wrapper for heat and salt transport calculations to calculate them for each basin
487      !! Called from all advection and/or diffusion routines
488      !!----------------------------------------------------------------------
489      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
490      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
491      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
492      INTEGER                                        :: jn    !
493
494      IF( cptr == 'adv' ) THEN
495         IF( ktra == jp_tem )  htr_adv(:,1) = ptr_sj( pva(:,:,:) )
496         IF( ktra == jp_sal )  str_adv(:,1) = ptr_sj( pva(:,:,:) )
497      ENDIF
498      IF( cptr == 'ldf' ) THEN
499         IF( ktra == jp_tem )  htr_ldf(:,1) = ptr_sj( pva(:,:,:) )
500         IF( ktra == jp_sal )  str_ldf(:,1) = ptr_sj( pva(:,:,:) )
501      ENDIF
502      IF( cptr == 'eiv' ) THEN
503         IF( ktra == jp_tem )  htr_eiv(:,1) = ptr_sj( pva(:,:,:) )
504         IF( ktra == jp_sal )  str_eiv(:,1) = ptr_sj( pva(:,:,:) )
505      ENDIF
506      !
507      IF( ln_subbas ) THEN
508         !
509         IF( cptr == 'adv' ) THEN
510             IF( ktra == jp_tem ) THEN
511                DO jn = 2, nptr
512                   htr_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
513                END DO
514             ENDIF
515             IF( ktra == jp_sal ) THEN
516                DO jn = 2, nptr
517                   str_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
518                END DO
519             ENDIF
520         ENDIF
521         IF( cptr == 'ldf' ) THEN
522             IF( ktra == jp_tem ) THEN
523                DO jn = 2, nptr
524                    htr_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
525                 END DO
526             ENDIF
527             IF( ktra == jp_sal ) THEN
528                DO jn = 2, nptr
529                   str_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
530                END DO
531             ENDIF
532         ENDIF
533         IF( cptr == 'eiv' ) THEN
534             IF( ktra == jp_tem ) THEN
535                DO jn = 2, nptr
536                    htr_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
537                 END DO
538             ENDIF
539             IF( ktra == jp_sal ) THEN
540                DO jn = 2, nptr
541                   str_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
542                END DO
543             ENDIF
544         ENDIF
545         !
546      ENDIF
547   END SUBROUTINE dia_ptr_hst
548
549
550   FUNCTION dia_ptr_alloc()
551      !!----------------------------------------------------------------------
552      !!                    ***  ROUTINE dia_ptr_alloc  ***
553      !!----------------------------------------------------------------------
554      INTEGER               ::   dia_ptr_alloc   ! return value
555      INTEGER, DIMENSION(3) ::   ierr
556      !!----------------------------------------------------------------------
557      ierr(:) = 0
558      !
559      ALLOCATE( btmsk(jpi,jpj,nptr) ,              &
560         &      htr_adv(jpj,nptr) , str_adv(jpj,nptr) ,   &
561         &      htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) ,   &
562         &      htr_ove(jpj,nptr) , str_ove(jpj,nptr) ,   &
563         &      htr_btr(jpj,nptr) , str_btr(jpj,nptr) ,   &
564         &      htr_ldf(jpj,nptr) , str_ldf(jpj,nptr) , STAT=ierr(1)  )
565         !
566      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2))
567      !
568      ALLOCATE( btm30(jpi,jpj), STAT=ierr(3)  )
569
570         !
571      dia_ptr_alloc = MAXVAL( ierr )
572      IF(lk_mpp)   CALL mpp_sum( dia_ptr_alloc )
573      !
574   END FUNCTION dia_ptr_alloc
575
576
577   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval )
578      !!----------------------------------------------------------------------
579      !!                    ***  ROUTINE ptr_sj_3d  ***
580      !!
581      !! ** Purpose :   i-k sum computation of a j-flux array
582      !!
583      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
584      !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
585      !!
586      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
587      !!----------------------------------------------------------------------
588      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pva   ! mask flux array at V-point
589      REAL(wp), INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
590      !
591      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
592      INTEGER                  ::   ijpj         ! ???
593      REAL(wp), POINTER, DIMENSION(:) :: p_fval  ! function value
594      !!--------------------------------------------------------------------
595      !
596      p_fval => p_fval1d
597
598      ijpj = jpj
599      p_fval(:) = 0._wp
600      IF( PRESENT( pmsk ) ) THEN
601         DO jk = 1, jpkm1
602            DO jj = 2, jpjm1
603               DO ji = fs_2, fs_jpim1   ! Vector opt.
604                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) * pmsk(ji,jj)
605               END DO
606            END DO
607         END DO
608      ELSE
609         DO jk = 1, jpkm1
610            DO jj = 2, jpjm1
611               DO ji = fs_2, fs_jpim1   ! Vector opt.
612                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) 
613               END DO
614            END DO
615         END DO
616      ENDIF
617#if defined key_mpp_mpi
618      IF(lk_mpp)   CALL mpp_sum( p_fval, ijpj, ncomm_znl)
619#endif
620      !
621   END FUNCTION ptr_sj_3d
622
623
624   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval )
625      !!----------------------------------------------------------------------
626      !!                    ***  ROUTINE ptr_sj_2d  ***
627      !!
628      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array
629      !!
630      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
631      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
632      !!
633      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
634      !!----------------------------------------------------------------------
635      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pva   ! mask flux array at V-point
636      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
637      !
638      INTEGER                  ::   ji,jj       ! dummy loop arguments
639      INTEGER                  ::   ijpj        ! ???
640      REAL(wp), POINTER, DIMENSION(:) :: p_fval ! function value
641      !!--------------------------------------------------------------------
642      !
643      p_fval => p_fval1d
644
645      ijpj = jpj
646      p_fval(:) = 0._wp
647      IF( PRESENT( pmsk ) ) THEN
648         DO jj = 2, jpjm1
649            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
650               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) * pmsk(ji,jj)
651            END DO
652         END DO
653      ELSE
654         DO jj = 2, jpjm1
655            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
656               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj)
657            END DO
658         END DO
659      ENDIF
660#if defined key_mpp_mpi
661      CALL mpp_sum( p_fval, ijpj, ncomm_znl )
662#endif
663      !
664   END FUNCTION ptr_sj_2d
665
666
667   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
668      !!----------------------------------------------------------------------
669      !!                    ***  ROUTINE ptr_sjk  ***
670      !!
671      !! ** Purpose :   i-sum computation of an array
672      !!
673      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
674      !!
675      !! ** Action  : - p_fval: i-mean poleward flux of pva
676      !!----------------------------------------------------------------------
677      !!
678      IMPLICIT none
679      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pta    ! mask flux array at V-point
680      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask
681      !!
682      INTEGER                           :: ji, jj, jk ! dummy loop arguments
683      REAL(wp), POINTER, DIMENSION(:,:) :: p_fval     ! return function value
684#if defined key_mpp_mpi
685      INTEGER, DIMENSION(1) ::   ish
686      INTEGER, DIMENSION(2) ::   ish2
687      INTEGER               ::   ijpjjpk
688      REAL(wp), DIMENSION(jpj*jpk) ::   zwork    ! mask flux array at V-point
689#endif
690      !!--------------------------------------------------------------------
691      !
692      p_fval => p_fval2d
693
694      p_fval(:,:) = 0._wp
695      !
696      IF( PRESENT( pmsk ) ) THEN
697         DO jk = 1, jpkm1
698            DO jj = 2, jpjm1
699!!gm here, use of tmask_i  ==> no need of loop over nldi, nlei....
700               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
701                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj)
702               END DO
703            END DO
704         END DO
705      ELSE
706         DO jk = 1, jpkm1
707            DO jj = 2, jpjm1
708               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
709                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * tmask_i(ji,jj)
710               END DO
711            END DO
712         END DO
713      END IF
714      !
715#if defined key_mpp_mpi
716      ijpjjpk = jpj*jpk
717      ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk
718      zwork(1:ijpjjpk) = RESHAPE( p_fval, ish )
719      CALL mpp_sum( zwork, ijpjjpk, ncomm_znl )
720      p_fval(:,:) = RESHAPE( zwork, ish2 )
721#endif
722      !
723   END FUNCTION ptr_sjk
724
725
726   !!======================================================================
727END MODULE diaptr
Note: See TracBrowser for help on using the repository browser.