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/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diaptr.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago

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