source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DIA/diaptr.F90 @ 13519

Last change on this file since 13519 was 13519, checked in by hadcv, 7 months ago

Tiling for DIA routines called by TRA modules

  • Property svn:keywords set to Id
File size: 33.2 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   !!            4.0  ! 2010-08  ( C. Ethe, J. Deshayes ) Improvment
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   dia_ptr      : Poleward Transport Diagnostics module
17   !!   dia_ptr_init : Initialization, namelist read
18   !!   ptr_sjk      : "zonal" mean computation of a field - tracer or flux array
19   !!   ptr_sj       : "zonal" and vertical sum computation of a "meridional" flux array
20   !!                   (Generic interface to ptr_sj_3d, ptr_sj_2d)
21   !!----------------------------------------------------------------------
22   USE oce              ! ocean dynamics and active tracers
23   USE dom_oce          ! ocean space and time domain
24   ! TEMP: Possibly not necessary if using XIOS (if cumulative axis operations are possible)
25   USE domain, ONLY : dom_tile
26   USE phycst           ! physical constants
27   !
28   USE iom              ! IOM library
29   USE in_out_manager   ! I/O manager
30   USE lib_mpp          ! MPP library
31   USE timing           ! preformance summary
32
33   IMPLICIT NONE
34   PRIVATE
35
36   INTERFACE ptr_sum
37      MODULE PROCEDURE ptr_sum_3d, ptr_sum_2d
38   END INTERFACE
39
40   INTERFACE ptr_sj
41      MODULE PROCEDURE ptr_sj_3d, ptr_sj_2d
42   END INTERFACE
43
44   PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines
45   PUBLIC   ptr_sjk        !
46   PUBLIC   dia_ptr_init   ! call in memogcm
47   PUBLIC   dia_ptr        ! call in step module
48   PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines
49
50   !                                  !!** namelist  namptr  **
51   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hstr_adv, hstr_ldf, hstr_eiv   !: Heat/Salt TRansports(adv, diff, Bolus.)
52   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional)
53   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   pvtr_int, pzon_int             !: Other zonal integrals
54
55   LOGICAL , PUBLIC ::   l_diaptr        !: tracers  trend flag (set from namelist in trdini)
56   INTEGER, PARAMETER, PUBLIC ::   nptr = 5  ! (glo, atl, pac, ind, ipc)
57   INTEGER, PARAMETER         ::   jp_msk = 3
58   INTEGER, PARAMETER         ::   jp_vtr = 4
59
60   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
61   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rho0 x Cp)
62   REAL(wp) ::   rc_ggram = 1.e-9_wp   ! conversion from g    to Gg  (further x rho0)
63
64   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
65   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk34 ! mask out Southern Ocean (=0 south of 34°S)
66
67   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini)
68   
69   !! * Substitutions
70#  include "do_loop_substitute.h90"
71#  include "domzgr_substitute.h90"
72   !!----------------------------------------------------------------------
73   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
74   !! $Id$
75   !! Software governed by the CeCILL license (see ./LICENSE)
76   !!----------------------------------------------------------------------
77CONTAINS
78
79   ! TEMP: Most changes and some code in this module not necessary if using XIOS (subdomain support, axis operations)
80   SUBROUTINE dia_ptr( kt, Kmm, pvtr )
81      !!----------------------------------------------------------------------
82      !!                  ***  ROUTINE dia_ptr  ***
83      !!----------------------------------------------------------------------
84      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index     
85      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
86      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk)    , INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
87      !!----------------------------------------------------------------------
88      !
89      IF( ln_timing )   CALL timing_start('dia_ptr')
90
91      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init
92      !
93      IF( l_diaptr ) THEN
94         ! Calculate zonal integrals
95         IF( PRESENT( pvtr ) ) THEN
96            CALL dia_ptr_zint( Kmm, pvtr )
97         ELSE
98            CALL dia_ptr_zint( Kmm )
99         ENDIF
100
101         ! Calculate diagnostics only when zonal integrals have finished
102         IF( ntile == 0 .OR. ntile == nijtile ) CALL dia_ptr_iom(kt, Kmm, pvtr)
103      ENDIF
104
105      IF( ln_timing )   CALL timing_stop('dia_ptr')
106      !
107   END SUBROUTINE dia_ptr
108
109
110   SUBROUTINE dia_ptr_iom( kt, Kmm, pvtr )
111      !!----------------------------------------------------------------------
112      !!                  ***  ROUTINE dia_ptr_iom  ***
113      !!----------------------------------------------------------------------
114      !! ** Purpose : Calculate diagnostics and send to XIOS
115      !!----------------------------------------------------------------------
116      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index
117      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
118      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk)    , INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
119      !
120      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
121      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
122      REAL(wp), DIMENSION(jpj)      ::  zvsum, ztsum, zssum   ! 1D workspace
123      !
124      !overturning calculation
125      REAL(wp), DIMENSION(jpj,jpk,nptr) :: sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse
126      REAL(wp), DIMENSION(jpj,jpk,nptr) :: zt_jk, zs_jk ! i-mean T and S, j-Stream-Function
127
128      REAL(wp), DIMENSION(jpi,jpj,jpk,nptr)  :: z4d1, z4d2
129      REAL(wp), DIMENSION(jpi,jpj,nptr)      :: z3dtr ! i-mean T and S, j-Stream-Function
130      !!----------------------------------------------------------------------
131
132      IF( PRESENT( pvtr ) ) THEN
133         IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF
134            DO jn = 1, nptr                                    ! by sub-basins
135               z4d1(1,:,:,jn) =  pvtr_int(:,:,jp_vtr,jn)                  ! zonal cumulative effective transport excluding closed seas
136               DO jk = jpkm1, 1, -1
137                  z4d1(1,:,jk,jn) = z4d1(1,:,jk+1,jn) - z4d1(1,:,jk,jn)    ! effective j-Stream-Function (MSF)
138               END DO
139               DO ji = 1, jpi
140                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
141               ENDDO
142            END DO
143            CALL iom_put( 'zomsf', z4d1 * rc_sv )
144         ENDIF
145         IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
146            DO jn = 1, nptr
147               sjk(:,:,jn) = pvtr_int(:,:,jp_msk,jn)
148               r1_sjk(:,:,jn) = 0._wp
149               WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
150               ! i-mean T and S, j-Stream-Function, basin
151               zt_jk(:,:,jn) = pvtr_int(:,:,jp_tem,jn) * r1_sjk(:,:,jn)
152               zs_jk(:,:,jn) = pvtr_int(:,:,jp_sal,jn) * r1_sjk(:,:,jn)
153               v_msf(:,:,jn) = pvtr_int(:,:,jp_vtr,jn)
154               hstr_ove(:,jp_tem,jn) = SUM( v_msf(:,:,jn)*zt_jk(:,:,jn), 2 )
155               hstr_ove(:,jp_sal,jn) = SUM( v_msf(:,:,jn)*zs_jk(:,:,jn), 2 )
156               !
157            ENDDO
158            DO jn = 1, nptr
159               z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
160               DO ji = 1, jpi
161                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
162               ENDDO
163            ENDDO
164            CALL iom_put( 'sophtove', z3dtr )
165            DO jn = 1, nptr
166               z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
167               DO ji = 1, jpi
168                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
169               ENDDO
170            ENDDO
171            CALL iom_put( 'sopstove', z3dtr )
172         ENDIF
173
174         IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
175            ! Calculate barotropic heat and salt transport here
176            DO jn = 1, nptr
177               sjk(:,1,jn) = SUM( pvtr_int(:,:,jp_msk,jn), 2 )
178               r1_sjk(:,1,jn) = 0._wp
179               WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
180               !
181               zvsum(:) =    SUM( pvtr_int(:,:,jp_vtr,jn), 2 )
182               ztsum(:) =    SUM( pvtr_int(:,:,jp_tem,jn), 2 )
183               zssum(:) =    SUM( pvtr_int(:,:,jp_sal,jn), 2 )
184               hstr_btr(:,jp_tem,jn) = zvsum(:) * ztsum(:) * r1_sjk(:,1,jn)
185               hstr_btr(:,jp_sal,jn) = zvsum(:) * zssum(:) * r1_sjk(:,1,jn)
186               !
187            ENDDO
188            DO jn = 1, nptr
189               z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
190               DO ji = 1, jpi
191                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
192               ENDDO
193            ENDDO
194            CALL iom_put( 'sophtbtr', z3dtr )
195            DO jn = 1, nptr
196               z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
197               DO ji = 1, jpi
198                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
199               ENDDO
200            ENDDO
201            CALL iom_put( 'sopstbtr', z3dtr )
202         ENDIF 
203         !
204         hstr_ove(:,:,:) = 0._wp       ! Zero before next timestep
205         hstr_btr(:,:,:) = 0._wp
206         pvtr_int(:,:,:,:) = 0._wp
207      ELSE
208         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface
209            !
210            DO jn = 1, nptr
211               z4d1(1,:,:,jn) = pzon_int(:,:,jp_msk,jn)
212               DO ji = 2, jpi
213                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
214               ENDDO
215            ENDDO
216            CALL iom_put( 'zosrf', z4d1 )
217            !
218            DO jn = 1, nptr
219               z4d2(1,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
220               DO ji = 2, jpi
221                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
222               ENDDO
223            ENDDO
224            CALL iom_put( 'zotem', z4d2 )
225            !
226            DO jn = 1, nptr
227               z4d2(1,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
228               DO ji = 2, jpi
229                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
230               ENDDO
231            ENDDO
232            CALL iom_put( 'zosal', z4d2 )
233            !
234         ENDIF
235         !
236         !                                ! Advective and diffusive heat and salt transport
237         IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN 
238            !
239            DO jn = 1, nptr
240               z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
241               DO ji = 1, jpi
242                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
243               ENDDO
244            ENDDO
245            CALL iom_put( 'sophtadv', z3dtr )
246            DO jn = 1, nptr
247               z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
248               DO ji = 1, jpi
249                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
250               ENDDO
251            ENDDO
252            CALL iom_put( 'sopstadv', z3dtr )
253         ENDIF
254         !
255         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN 
256            !
257            DO jn = 1, nptr
258               z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
259               DO ji = 1, jpi
260                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
261               ENDDO
262            ENDDO
263            CALL iom_put( 'sophtldf', z3dtr )
264            DO jn = 1, nptr
265               z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
266               DO ji = 1, jpi
267                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
268               ENDDO
269            ENDDO
270            CALL iom_put( 'sopstldf', z3dtr )
271         ENDIF
272         !
273         IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN 
274            !
275            DO jn = 1, nptr
276               z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
277               DO ji = 1, jpi
278                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
279               ENDDO
280            ENDDO
281            CALL iom_put( 'sophteiv', z3dtr )
282            DO jn = 1, nptr
283               z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
284               DO ji = 1, jpi
285                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
286               ENDDO
287            ENDDO
288            CALL iom_put( 'sopsteiv', z3dtr )
289         ENDIF
290         !
291         IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
292             DO jn = 1, nptr
293                z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
294                DO ji = 1, jpi
295                   z3dtr(ji,:,jn) = z3dtr(1,:,jn)
296                ENDDO
297             ENDDO
298             CALL iom_put( 'sophtvtr', z3dtr )
299             DO jn = 1, nptr
300               z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
301               DO ji = 1, jpi
302                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
303               ENDDO
304            ENDDO
305            CALL iom_put( 'sopstvtr', z3dtr )
306         ENDIF
307         !
308         ! TEMP: Possibly not necessary if using XIOS (if cumulative axis operations are possible)
309         ! TODO: NOT TESTED- hangs on iom_get_var
310         IF( iom_use( 'uocetr_vsum_cumul' ) ) THEN
311            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
312            CALL iom_get_var(  'uocetr_vsum_op', z2d ) ! get uocetr_vsum_op from xml
313            z2d(:,:) = ptr_ci_2d( z2d(:,:) ) 
314            CALL iom_put( 'uocetr_vsum_cumul', z2d )
315            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )   ! Revert to tile domain
316         ENDIF
317         !
318         hstr_adv(:,:,:) = 0._wp       ! Zero before next timestep
319         hstr_ldf(:,:,:) = 0._wp
320         hstr_eiv(:,:,:) = 0._wp
321         hstr_vtr(:,:,:) = 0._wp
322         pzon_int(:,:,:,:) = 0._wp
323      ENDIF
324   END SUBROUTINE dia_ptr_iom
325
326
327   SUBROUTINE dia_ptr_zint( Kmm, pvtr )
328      !!----------------------------------------------------------------------
329      !!                    ***  ROUTINE dia_ptr_zint ***
330      !!----------------------------------------------------------------------
331      !! ** Purpose : i and i-k sum operations on arrays
332      !!
333      !! ** Method  : - Call ptr_sjk (i sum) or ptr_sj (i-k sum) to perform the sum operation
334      !!              - Call ptr_sum to add this result to the sum over tiles
335      !!
336      !! ** Action  : pvtr_int - terms for volume streamfunction, heat/salt transport barotropic/overturning terms
337      !!              pzon_int - terms for i mean temperature/salinity
338      !!----------------------------------------------------------------------
339      INTEGER                     , INTENT(in)           :: Kmm          ! time level index
340      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk), INTENT(in), OPTIONAL :: pvtr         ! j-effective transport
341      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk)                       :: zmask        ! 3D workspace
342      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk,jpts)                  :: zts          ! 4D workspace
343      REAL(wp), DIMENSION(ST_1Dj(nn_hls),jpk,nptr)                 :: sjk, v_msf   ! Zonal sum: i-k surface area, j-effective transport
344      REAL(wp), DIMENSION(ST_1Dj(nn_hls),jpk,nptr)                 :: zt_jk, zs_jk ! Zonal sum: i-k surface area * (T, S)
345      REAL(wp)                                           :: zsfc, zvfc   ! i-k surface area
346      INTEGER  ::   ji, jj, jk, jn                                       ! dummy loop indices
347      !!----------------------------------------------------------------------
348
349      IF( PRESENT( pvtr ) ) THEN
350         ! i sum of effective j transport excluding closed seas
351         IF( iom_use( 'zomsf' ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
352            DO jn = 1, nptr
353               v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )
354            ENDDO
355
356            CALL ptr_sum( pvtr_int(:,:,jp_vtr,:), v_msf(:,:,:) )
357         ENDIF
358
359         ! i sum of j surface area, j surface area - temperature/salinity product on V grid
360         IF(  iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.   &
361            & iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
362            zmask(:,:,:) = 0._wp
363            zts(:,:,:,:) = 0._wp
364
365            DO_3D( 1, 0, 1, 1, 1, jpkm1 )
366               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
367               zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
368               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc !Tracers averaged onto V grid
369               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
370            END_3D
371
372            DO jn = 1, nptr
373               sjk(:,:,jn)   = ptr_sjk( zmask(:,:,:)     , btmsk(:,:,jn) )
374               zt_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
375               zs_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
376            ENDDO
377
378            CALL ptr_sum( pvtr_int(:,:,jp_msk,:), sjk(:,:,:)   )
379            CALL ptr_sum( pvtr_int(:,:,jp_tem,:), zt_jk(:,:,:) )
380            CALL ptr_sum( pvtr_int(:,:,jp_sal,:), zs_jk(:,:,:) )
381         ENDIF
382      ELSE
383         ! i sum of j surface area - temperature/salinity product on T grid
384         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN
385            zmask(:,:,:) = 0._wp
386            zts(:,:,:,:) = 0._wp
387
388            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
389               zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm)
390               zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
391               zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc
392               zts(ji,jj,jk,jp_sal) = ts(ji,jj,jk,jp_sal,Kmm) * zsfc
393            END_3D
394
395            DO jn = 1, nptr
396               sjk(:,:,jn)   = ptr_sjk( zmask(:,:,:)     , btmsk(:,:,jn) )
397               zt_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
398               zs_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
399            ENDDO
400
401            CALL ptr_sum( pzon_int(:,:,jp_msk,:), sjk(:,:,:)   )
402            CALL ptr_sum( pzon_int(:,:,jp_tem,:), zt_jk(:,:,:) )
403            CALL ptr_sum( pzon_int(:,:,jp_sal,:), zs_jk(:,:,:) )
404         ENDIF
405
406         ! i-k sum of j surface area - temperature/salinity product on V grid
407         IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
408            zts(:,:,:,:) = 0._wp
409
410            DO_3D( 1, 0, 1, 1, 1, jpkm1 )
411               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
412               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid
413               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
414            END_3D
415
416            CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) )
417            CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) )
418         ENDIF
419      ENDIF
420   END SUBROUTINE dia_ptr_zint
421
422
423   SUBROUTINE dia_ptr_init
424      !!----------------------------------------------------------------------
425      !!                  ***  ROUTINE dia_ptr_init  ***
426      !!                   
427      !! ** Purpose :   Initialization, namelist read
428      !!----------------------------------------------------------------------
429      INTEGER ::  inum, jn           ! local integers
430      !!
431      REAL(wp), DIMENSION(jpi,jpj) :: zmsk
432      !!----------------------------------------------------------------------
433
434      l_diaptr = .FALSE.
435      IF(   iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  &
436         &  iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  &
437         &  iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  &
438         &  iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  & 
439         &  iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  &
440         &  iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) )  l_diaptr  = .TRUE.
441
442 
443      IF(lwp) THEN                     ! Control print
444         WRITE(numout,*)
445         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
446         WRITE(numout,*) '~~~~~~~~~~~~'
447         WRITE(numout,*) '   Namelist namptr : set ptr parameters'
448         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr
449      ENDIF
450
451      IF( l_diaptr ) THEN 
452         !
453         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
454
455         rc_pwatt = rc_pwatt * rho0_rcp          ! conversion from K.s-1 to PetaWatt
456         rc_ggram = rc_ggram * rho0              ! conversion from m3/s to Gg/s
457
458         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
459
460         btmsk(:,:,:) = 0._wp
461         btmsk(:,:,1) = tmask_i(:,:)
462         CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  )
463         CALL iom_get( inum, jpdom_global, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
464         CALL iom_get( inum, jpdom_global, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
465         CALL iom_get( inum, jpdom_global, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
466         CALL iom_close( inum )
467         btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
468         DO jn = 2, nptr
469            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
470         END DO
471         ! JD : modification so that overturning streamfunction is available in Atlantic at 34S to compare with observations
472         WHERE( gphit(:,:)*tmask_i(:,:) < -34._wp)
473           zmsk(:,:) = 0._wp      ! mask out Southern Ocean
474         ELSE WHERE                 
475           zmsk(:,:) = ssmask(:,:)
476         END WHERE
477         btmsk34(:,:,1) = btmsk(:,:,1)                 
478         DO jn = 2, nptr
479            btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:)               ! interior domain only
480         ENDDO
481
482         ! Initialise arrays to zero because diatpr is called before they are first calculated
483         ! Note that this means diagnostics will not be exactly correct when model run is restarted.
484         hstr_adv(:,:,:) = 0._wp           
485         hstr_ldf(:,:,:) = 0._wp           
486         hstr_eiv(:,:,:) = 0._wp           
487         hstr_ove(:,:,:) = 0._wp           
488         hstr_btr(:,:,:) = 0._wp           !
489         hstr_vtr(:,:,:) = 0._wp           !
490         pvtr_int(:,:,:,:) = 0._wp
491         pzon_int(:,:,:,:) = 0._wp
492         !
493         ll_init = .FALSE.
494         !
495      ENDIF 
496      !
497   END SUBROUTINE dia_ptr_init
498
499
500   SUBROUTINE dia_ptr_hst( ktra, cptr, pvflx ) 
501      !!----------------------------------------------------------------------
502      !!                    ***  ROUTINE dia_ptr_hst ***
503      !!----------------------------------------------------------------------
504      !! Wrapper for heat and salt transport calculations to calculate them for each basin
505      !! Called from all advection and/or diffusion routines
506      !!----------------------------------------------------------------------
507      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
508      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
509      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk)    , INTENT(in)   :: pvflx ! 3D input array of advection/diffusion
510      REAL(wp), DIMENSION(ST_1Dj(nn_hls),nptr)                 :: zsj   !
511      INTEGER                                        :: jn    !
512
513      DO jn = 1, nptr
514         zsj(:,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
515      ENDDO
516      !
517      IF( cptr == 'adv' ) THEN
518         IF( ktra == jp_tem )  CALL ptr_sum( hstr_adv(:,jp_tem,:), zsj(:,:) )
519         IF( ktra == jp_sal )  CALL ptr_sum( hstr_adv(:,jp_sal,:), zsj(:,:) )
520      ELSE IF( cptr == 'ldf' ) THEN
521         IF( ktra == jp_tem )  CALL ptr_sum( hstr_ldf(:,jp_tem,:), zsj(:,:) )
522         IF( ktra == jp_sal )  CALL ptr_sum( hstr_ldf(:,jp_sal,:), zsj(:,:) )
523      ELSE IF( cptr == 'eiv' ) THEN
524         IF( ktra == jp_tem )  CALL ptr_sum( hstr_eiv(:,jp_tem,:), zsj(:,:) )
525         IF( ktra == jp_sal )  CALL ptr_sum( hstr_eiv(:,jp_sal,:), zsj(:,:) )
526      ELSE IF( cptr == 'vtr' ) THEN
527         IF( ktra == jp_tem )  CALL ptr_sum( hstr_vtr(:,jp_tem,:), zsj(:,:) )
528         IF( ktra == jp_sal )  CALL ptr_sum( hstr_vtr(:,jp_sal,:), zsj(:,:) )
529      ENDIF
530      !
531   END SUBROUTINE dia_ptr_hst
532
533
534   SUBROUTINE ptr_sum_2d( phstr, pva )
535      !!----------------------------------------------------------------------
536      !!                    ***  ROUTINE ptr_sum_2d ***
537      !!----------------------------------------------------------------------
538      !! ** Purpose : Add two 2D arrays with (j,nptr) dimensions
539      !!
540      !! ** Method  : - phstr = phstr + pva
541      !!              - Call mpp_sum if the final tile
542      !!
543      !! ** Action  : phstr
544      !!----------------------------------------------------------------------
545      REAL(wp), DIMENSION(jpj,nptr) , INTENT(inout)         ::  phstr  !
546      REAL(wp), DIMENSION(ST_1Dj(nn_hls),nptr), INTENT(in)            ::  pva    !
547      INTEGER                                               ::  jj
548#if defined key_mpp_mpi
549      INTEGER, DIMENSION(1)          ::  ish1d
550      INTEGER, DIMENSION(2)          ::  ish2d
551      REAL(wp), DIMENSION(jpj*nptr)  ::  zwork
552#endif
553
554      DO jj = ntsj, ntej
555         phstr(jj,:) = phstr(jj,:)  + pva(jj,:)
556      END DO
557
558#if defined key_mpp_mpi
559      IF( ntile == 0 .OR. ntile == nijtile ) THEN
560         ish1d(1) = jpj*nptr
561         ish2d(1) = jpj ; ish2d(2) = nptr
562         zwork(:) = RESHAPE( phstr(:,:), ish1d )
563         CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
564         phstr(:,:) = RESHAPE( zwork, ish2d )
565      ENDIF
566#endif
567   END SUBROUTINE ptr_sum_2d
568
569
570   SUBROUTINE ptr_sum_3d( phstr, pva )
571      !!----------------------------------------------------------------------
572      !!                    ***  ROUTINE ptr_sum_3d ***
573      !!----------------------------------------------------------------------
574      !! ** Purpose : Add two 3D arrays with (j,k,nptr) dimensions
575      !!
576      !! ** Method  : - phstr = phstr + pva
577      !!              - Call mpp_sum if the final tile
578      !!
579      !! ** Action  : phstr
580      !!----------------------------------------------------------------------
581      REAL(wp), DIMENSION(jpj,jpk,nptr) , INTENT(inout)     ::  phstr  !
582      REAL(wp), DIMENSION(ST_1Dj(nn_hls),jpk,nptr), INTENT(in)        ::  pva    !
583      INTEGER                                               ::  jj, jk
584#if defined key_mpp_mpi
585      INTEGER, DIMENSION(1)              ::  ish1d
586      INTEGER, DIMENSION(3)              ::  ish3d
587      REAL(wp), DIMENSION(jpj*jpk*nptr)  ::  zwork
588#endif
589
590      DO jk = 1, jpk
591         DO jj = ntsj, ntej
592            phstr(jj,jk,:) = phstr(jj,jk,:)  + pva(jj,jk,:)
593         END DO
594      END DO
595
596#if defined key_mpp_mpi
597      IF( ntile == 0 .OR. ntile == nijtile ) THEN
598         ish1d(1) = jpj*jpk*nptr
599         ish3d(1) = jpj ; ish3d(2) = jpk ; ish3d(3) = nptr
600         zwork(:) = RESHAPE( phstr(:,:,:), ish1d )
601         CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
602         phstr(:,:,:) = RESHAPE( zwork, ish3d )
603      ENDIF
604#endif
605   END SUBROUTINE ptr_sum_3d
606
607
608   FUNCTION dia_ptr_alloc()
609      !!----------------------------------------------------------------------
610      !!                    ***  ROUTINE dia_ptr_alloc  ***
611      !!----------------------------------------------------------------------
612      INTEGER               ::   dia_ptr_alloc   ! return value
613      INTEGER, DIMENSION(2) ::   ierr
614      !!----------------------------------------------------------------------
615      ierr(:) = 0
616      !
617      IF( .NOT. ALLOCATED( btmsk ) ) THEN
618         ALLOCATE( btmsk(jpi,jpj,nptr)    , btmsk34(jpi,jpj,nptr),   &
619            &      hstr_adv(jpj,jpts,nptr), hstr_eiv(jpj,jpts,nptr), &
620            &      hstr_ove(jpj,jpts,nptr), hstr_btr(jpj,jpts,nptr), &
621            &      hstr_ldf(jpj,jpts,nptr), hstr_vtr(jpj,jpts,nptr), STAT=ierr(1)  )
622            !
623         ALLOCATE( pvtr_int(jpj,jpk,jpts+2,nptr), &
624            &      pzon_int(jpj,jpk,jpts+1,nptr), STAT=ierr(2) )
625         !
626         dia_ptr_alloc = MAXVAL( ierr )
627         CALL mpp_sum( 'diaptr', dia_ptr_alloc )
628      ENDIF
629      !
630   END FUNCTION dia_ptr_alloc
631
632
633   FUNCTION ptr_sj_3d( pvflx, pmsk )   RESULT ( p_fval )
634      !!----------------------------------------------------------------------
635      !!                    ***  ROUTINE ptr_sj_3d  ***
636      !!
637      !! ** Purpose :   i-k sum computation of a j-flux array
638      !!
639      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
640      !!              pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
641      !!
642      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
643      !!----------------------------------------------------------------------
644      REAL(wp), INTENT(in), DIMENSION(ST_2D(nn_hls),jpk)  ::   pvflx  ! mask flux array at V-point
645      REAL(wp), INTENT(in), DIMENSION(jpi,jpj)  ::   pmsk   ! Optional 2D basin mask
646      !
647      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
648      REAL(wp), DIMENSION(ST_1Dj(nn_hls)) :: p_fval  ! function value
649      !!--------------------------------------------------------------------
650      !
651      p_fval(:) = 0._wp
652      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
653         p_fval(jj) = p_fval(jj) + pvflx(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
654      END_3D
655   END FUNCTION ptr_sj_3d
656
657
658   FUNCTION ptr_sj_2d( pvflx, pmsk )   RESULT ( p_fval )
659      !!----------------------------------------------------------------------
660      !!                    ***  ROUTINE ptr_sj_2d  ***
661      !!
662      !! ** Purpose :   "zonal" and vertical sum computation of a j-flux array
663      !!
664      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
665      !!      pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
666      !!
667      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
668      !!----------------------------------------------------------------------
669      REAL(wp) , INTENT(in), DIMENSION(ST_2D(nn_hls))     ::   pvflx  ! mask flux array at V-point
670      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
671      !
672      INTEGER                  ::   ji,jj       ! dummy loop arguments
673      REAL(wp), DIMENSION(ST_1Dj(nn_hls)) :: p_fval ! function value
674      !!--------------------------------------------------------------------
675      !
676      p_fval(:) = 0._wp
677      DO_2D( 0, 0, 0, 0 )
678         p_fval(jj) = p_fval(jj) + pvflx(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj)
679      END_2D
680   END FUNCTION ptr_sj_2d
681
682   FUNCTION ptr_ci_2d( pva )   RESULT ( p_fval )
683      !!----------------------------------------------------------------------
684      !!                    ***  ROUTINE ptr_ci_2d  ***
685      !!
686      !! ** Purpose :   "meridional" cumulated sum computation of a j-flux array
687      !!
688      !! ** Method  : - j cumulated sum of pva using the interior 2D vmask (umask_i).
689      !!
690      !! ** Action  : - p_fval: j-cumulated sum of pva
691      !!----------------------------------------------------------------------
692      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)  ::   pva   ! mask flux array at V-point
693      !
694      INTEGER                  ::   ji,jj,jc       ! dummy loop arguments
695      INTEGER                  ::   ijpj        ! ???
696      REAL(wp), DIMENSION(jpi,jpj) :: p_fval ! function value
697      !!--------------------------------------------------------------------
698      !
699      ijpj = jpj  ! ???
700      p_fval(:,:) = 0._wp
701      DO jc = 1, jpnj ! looping over all processors in j axis
702         DO_2D( 0, 0, 0, 0 )
703            p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj)
704         END_2D
705         CALL lbc_lnk( 'diaptr', p_fval, 'U', -1.0_wp )
706      END DO
707      !
708   END FUNCTION ptr_ci_2d
709
710
711
712   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
713      !!----------------------------------------------------------------------
714      !!                    ***  ROUTINE ptr_sjk  ***
715      !!
716      !! ** Purpose :   i-sum computation of an array
717      !!
718      !! ** Method  : - i-sum of field using the interior 2D vmask (pmsk).
719      !!
720      !! ** Action  : - p_fval: i-sum of masked field
721      !!----------------------------------------------------------------------
722      !!
723      IMPLICIT none
724      REAL(wp) , INTENT(in), DIMENSION(ST_2D(nn_hls),jpk) ::   pta    ! mask flux array at V-point
725      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
726      !!
727      INTEGER                           :: ji, jj, jk ! dummy loop arguments
728      REAL(wp), DIMENSION(ST_1Dj(nn_hls),jpk) :: p_fval     ! return function value
729      !!--------------------------------------------------------------------
730      !
731      p_fval(:,:) = 0._wp
732      !
733      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
734         p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
735      END_3D
736   END FUNCTION ptr_sjk
737
738
739   !!======================================================================
740END MODULE diaptr
Note: See TracBrowser for help on using the repository browser.