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 NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DIA – NEMO

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

Last change on this file since 13741 was 13741, checked in by hadcv, 4 years ago

#2365: Merge in trunk changes to [13688] for src/cfgs

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