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/trunk/src/OCE/DIA – NEMO

source: NEMO/trunk/src/OCE/DIA/diaptr.F90 @ 15004

Last change on this file since 15004 was 14834, checked in by hadcv, 3 years ago

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