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 @ 14229

Last change on this file since 14229 was 14229, checked in by smasson, 3 years ago

trunk: replace key_mpp_mpi by the opposite key key_mpi_off
  • Property svn:keywords set to Id
File size: 34.0 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   SUBROUTINE dia_ptr( kt, Kmm, pvtr )
74      !!----------------------------------------------------------------------
75      !!                  ***  ROUTINE dia_ptr  ***
76      !!----------------------------------------------------------------------
77      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index
78      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
79      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
80      !!----------------------------------------------------------------------
81      !
82      IF( ln_timing )   CALL timing_start('dia_ptr')
83
84      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init    ! -> will define l_diaptr and nbasin
85      !
86      IF( l_diaptr ) THEN
87         ! Calculate zonal integrals
88         IF( PRESENT( pvtr ) ) THEN
89            CALL dia_ptr_zint( Kmm, pvtr )
90         ELSE
91            CALL dia_ptr_zint( Kmm )
92         ENDIF
93
94         ! Calculate diagnostics only when zonal integrals have finished
95         IF( ntile == 0 .OR. ntile == nijtile ) CALL dia_ptr_iom(kt, Kmm, pvtr)
96      ENDIF
97
98      IF( ln_timing )   CALL timing_stop('dia_ptr')
99      !
100   END SUBROUTINE dia_ptr
101
102
103   SUBROUTINE dia_ptr_iom( kt, Kmm, pvtr )
104      !!----------------------------------------------------------------------
105      !!                  ***  ROUTINE dia_ptr_iom  ***
106      !!----------------------------------------------------------------------
107      !! ** Purpose : Calculate diagnostics and send to XIOS
108      !!----------------------------------------------------------------------
109      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index
110      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
111      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
112      !
113      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
114      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
115      REAL(wp), DIMENSION(jpj)      ::  zvsum, ztsum, zssum   ! 1D workspace
116      !
117      !overturning calculation
118      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse
119      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   zt_jk, zs_jk        ! i-mean T and S, j-Stream-Function
120
121      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   z4d1, z4d2
122      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   z3dtr
123      !!----------------------------------------------------------------------
124      !
125      ALLOCATE( z3dtr(jpi,jpj,nbasin) )
126
127      IF( PRESENT( pvtr ) ) THEN
128         IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF
129            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin) )
130            !
131            DO jn = 1, nbasin                                    ! by sub-basins
132               z4d1(1,:,:,jn) =  pvtr_int(:,:,jp_vtr,jn)                  ! zonal cumulative effective transport excluding closed seas
133               DO jk = jpkm1, 1, -1
134                  z4d1(1,:,jk,jn) = z4d1(1,:,jk+1,jn) - z4d1(1,:,jk,jn)    ! effective j-Stream-Function (MSF)
135               END DO
136               DO ji = 2, jpi
137                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
138               ENDDO
139            END DO
140            CALL iom_put( 'zomsf', z4d1 * rc_sv )
141            !
142            DEALLOCATE( z4d1 )
143         ENDIF
144         IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
145            ALLOCATE( sjk(jpj,jpk,nbasin), r1_sjk(jpj,jpk,nbasin), v_msf(jpj,jpk,nbasin),   &
146               &      zt_jk(jpj,jpk,nbasin), zs_jk(jpj,jpk,nbasin) )
147            !
148            DO jn = 1, nbasin
149               sjk(:,:,jn) = pvtr_int(:,:,jp_msk,jn)
150               r1_sjk(:,:,jn) = 0._wp
151               WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
152               ! i-mean T and S, j-Stream-Function, basin
153               zt_jk(:,:,jn) = pvtr_int(:,:,jp_tem,jn) * r1_sjk(:,:,jn)
154               zs_jk(:,:,jn) = pvtr_int(:,:,jp_sal,jn) * r1_sjk(:,:,jn)
155               v_msf(:,:,jn) = pvtr_int(:,:,jp_vtr,jn)
156               hstr_ove(:,jp_tem,jn) = SUM( v_msf(:,:,jn)*zt_jk(:,:,jn), 2 )
157               hstr_ove(:,jp_sal,jn) = SUM( v_msf(:,:,jn)*zs_jk(:,:,jn), 2 )
158               !
159            ENDDO
160            DO jn = 1, nbasin
161               z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
162               DO ji = 2, jpi
163                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
164               ENDDO
165            ENDDO
166            CALL iom_put( 'sophtove', z3dtr )
167            DO jn = 1, nbasin
168               z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
169               DO ji = 2, jpi
170                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
171               ENDDO
172            ENDDO
173            CALL iom_put( 'sopstove', z3dtr )
174            !
175            DEALLOCATE( sjk, r1_sjk, v_msf, zt_jk, zs_jk )
176         ENDIF
177
178         IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
179            ! Calculate barotropic heat and salt transport here
180            ALLOCATE( sjk(jpj,1,nbasin), r1_sjk(jpj,1,nbasin) )
181            !
182            DO jn = 1, nbasin
183               sjk(:,1,jn) = SUM( pvtr_int(:,:,jp_msk,jn), 2 )
184               r1_sjk(:,1,jn) = 0._wp
185               WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
186               !
187               zvsum(:) =    SUM( pvtr_int(:,:,jp_vtr,jn), 2 )
188               ztsum(:) =    SUM( pvtr_int(:,:,jp_tem,jn), 2 )
189               zssum(:) =    SUM( pvtr_int(:,:,jp_sal,jn), 2 )
190               hstr_btr(:,jp_tem,jn) = zvsum(:) * ztsum(:) * r1_sjk(:,1,jn)
191               hstr_btr(:,jp_sal,jn) = zvsum(:) * zssum(:) * r1_sjk(:,1,jn)
192               !
193            ENDDO
194            DO jn = 1, nbasin
195               z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
196               DO ji = 2, jpi
197                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
198               ENDDO
199            ENDDO
200            CALL iom_put( 'sophtbtr', z3dtr )
201            DO jn = 1, nbasin
202               z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
203               DO ji = 2, jpi
204                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
205               ENDDO
206            ENDDO
207            CALL iom_put( 'sopstbtr', z3dtr )
208            !
209            DEALLOCATE( sjk, r1_sjk )
210         ENDIF
211         !
212         hstr_ove(:,:,:) = 0._wp       ! Zero before next timestep
213         hstr_btr(:,:,:) = 0._wp
214         pvtr_int(:,:,:,:) = 0._wp
215      ELSE
216         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface
217            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin), z4d2(jpi,jpj,jpk,nbasin) )
218            !
219            DO jn = 1, nbasin
220               z4d1(1,:,:,jn) = pzon_int(:,:,jp_msk,jn)
221               DO ji = 2, jpi
222                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
223               ENDDO
224            ENDDO
225            CALL iom_put( 'zosrf', z4d1 )
226            !
227            DO jn = 1, nbasin
228               z4d2(1,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
229               DO ji = 2, jpi
230                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
231               ENDDO
232            ENDDO
233            CALL iom_put( 'zotem', z4d2 )
234            !
235            DO jn = 1, nbasin
236               z4d2(1,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
237               DO ji = 2, jpi
238                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
239               ENDDO
240            ENDDO
241            CALL iom_put( 'zosal', z4d2 )
242            !
243            DEALLOCATE( z4d1, z4d2 )
244         ENDIF
245         !
246         !                                ! Advective and diffusive heat and salt transport
247         IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN
248            !
249            DO jn = 1, nbasin
250               z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
251               DO ji = 2, jpi
252                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
253               ENDDO
254            ENDDO
255            CALL iom_put( 'sophtadv', z3dtr )
256            DO jn = 1, nbasin
257               z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
258               DO ji = 2, jpi
259                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
260               ENDDO
261            ENDDO
262            CALL iom_put( 'sopstadv', z3dtr )
263         ENDIF
264         !
265         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN
266            !
267            DO jn = 1, nbasin
268               z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
269               DO ji = 2, jpi
270                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
271               ENDDO
272            ENDDO
273            CALL iom_put( 'sophtldf', z3dtr )
274            DO jn = 1, nbasin
275               z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
276               DO ji = 2, jpi
277                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
278               ENDDO
279            ENDDO
280            CALL iom_put( 'sopstldf', z3dtr )
281         ENDIF
282         !
283         IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN
284            !
285            DO jn = 1, nbasin
286               z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
287               DO ji = 2, jpi
288                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
289               ENDDO
290            ENDDO
291            CALL iom_put( 'sophteiv', z3dtr )
292            DO jn = 1, nbasin
293               z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
294               DO ji = 2, jpi
295                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
296               ENDDO
297            ENDDO
298            CALL iom_put( 'sopsteiv', z3dtr )
299         ENDIF
300         !
301         IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
302             DO jn = 1, nbasin
303                z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
304                DO ji = 2, jpi
305                   z3dtr(ji,:,jn) = z3dtr(1,:,jn)
306                ENDDO
307             ENDDO
308             CALL iom_put( 'sophtvtr', z3dtr )
309             DO jn = 1, nbasin
310               z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
311               DO ji = 2, jpi
312                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
313               ENDDO
314            ENDDO
315            CALL iom_put( 'sopstvtr', z3dtr )
316         ENDIF
317         !
318         IF( iom_use( 'uocetr_vsum_cumul' ) ) THEN
319            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
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            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile )   ! Revert to tile domain
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(wp), DIMENSION(:,:,:), ALLOCATABLE            :: zmask        ! 3D workspace
353      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          :: zts          ! 4D workspace
354      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: sjk, v_msf   ! Zonal sum: i-k surface area, j-effective transport
355      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: zt_jk, zs_jk ! Zonal sum: i-k surface area * (T, S)
356      REAL(wp)                                           :: 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( 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(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: pvflx ! 3D input array of advection/diffusion
542      REAL(wp), 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(wp), DIMENSION(jpj,nbasin) , INTENT(inout)         ::  phstr  !
578      REAL(wp), 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(wp), 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( ntile == 0 .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(wp), DIMENSION(jpj,jpk,nbasin) , INTENT(inout)     ::  phstr  !
614      REAL(wp), 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(wp), 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( ntile == 0 .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(wp), INTENT(in), DIMENSION(A2D(nn_hls),jpk)  ::   pvflx  ! mask flux array at V-point
679      REAL(wp), INTENT(in), DIMENSION(jpi,jpj)  ::   pmsk   ! Optional 2D basin mask
680      !
681      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
682      REAL(wp), 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(wp), 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(wp) , 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(wp), 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(wp) , INTENT(in), DIMENSION(A2D(nn_hls),jpk) ::   pta    ! mask flux array at V-point
758      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
759      !!
760      INTEGER                           :: ji, jj, jk ! dummy loop arguments
761      REAL(wp), 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.