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/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/OCE/DIA/diaptr.F90 @ 13054

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

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