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

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DIA/diaptr.F90 @ 15513

Last change on this file since 15513 was 15513, checked in by techene, 3 years ago

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