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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diaptr.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 26.6 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 phycst           ! physical constants
25   !
26   USE iom              ! IOM library
27   USE in_out_manager   ! I/O manager
28   USE lib_mpp          ! MPP library
29   USE timing           ! preformance summary
30
31   IMPLICIT NONE
32   PRIVATE
33
34   INTERFACE ptr_sj
35      MODULE PROCEDURE ptr_sj_3d, ptr_sj_2d
36   END INTERFACE
37
38   PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines
39   PUBLIC   ptr_sjk        !
40   PUBLIC   dia_ptr_init   ! call in memogcm
41   PUBLIC   dia_ptr        ! call in step module
42   PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines
43
44   !                                  !!** namelist  namptr  **
45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_adv, hstr_ldf, hstr_eiv   !: Heat/Salt TRansports(adv, diff, Bolus.)
46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional)
47
48   LOGICAL , PUBLIC ::   l_diaptr        !: tracers  trend flag (set from namelist in trdini)
49   INTEGER, PARAMETER, PUBLIC ::   nptr = 5  ! (glo, atl, pac, ind, ipc)
50
51   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
52   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rau0 x Cp)
53   REAL(wp) ::   rc_ggram = 1.e-9_wp   ! conversion from g    to Gg  (further x rau0)
54
55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
56   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk34 ! mask out Southern Ocean (=0 south of 34°S)
57
58   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)   :: p_fval1d
59   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:) :: p_fval2d
60
61   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini)
62   !! * Substitutions
63#  include "vectopt_loop_substitute.h90"
64#  include "do_loop_substitute.h90"
65   !!----------------------------------------------------------------------
66   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
67   !! $Id$
68   !! Software governed by the CeCILL license (see ./LICENSE)
69   !!----------------------------------------------------------------------
70CONTAINS
71
72   SUBROUTINE dia_ptr( kt, Kmm, pvtr )
73      !!----------------------------------------------------------------------
74      !!                  ***  ROUTINE dia_ptr  ***
75      !!----------------------------------------------------------------------
76      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index     
77      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
78      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
79      !
80      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
81      REAL(wp) ::   zsfc,zvfc               ! local scalar
82      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
83      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace
84      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d    ! 3D workspace
85      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace
86      REAL(wp), DIMENSION(jpj)      ::  zvsum, ztsum, zssum   ! 1D workspace
87      !
88      !overturning calculation
89      REAL(wp), DIMENSION(jpj,jpk,nptr) :: sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse
90      REAL(wp), DIMENSION(jpj,jpk,nptr) :: zt_jk, zs_jk ! i-mean T and S, j-Stream-Function
91
92      REAL(wp), DIMENSION(jpi,jpj,jpk,nptr)  :: z4d1, z4d2
93      REAL(wp), DIMENSION(jpi,jpj,nptr)      :: z3dtr ! i-mean T and S, j-Stream-Function
94      !!----------------------------------------------------------------------
95      !
96      IF( ln_timing )   CALL timing_start('dia_ptr')
97
98      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init
99      !
100      IF( .NOT. l_diaptr )   RETURN
101
102      IF( PRESENT( pvtr ) ) THEN
103         IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF
104            DO jn = 1, nptr                                    ! by sub-basins
105               z4d1(1,:,:,jn) =  ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )  ! zonal cumulative effective transport excluding closed seas
106               DO jk = jpkm1, 1, -1 
107                  z4d1(1,:,jk,jn) = z4d1(1,:,jk+1,jn) - z4d1(1,:,jk,jn)    ! effective j-Stream-Function (MSF)
108               END DO
109               DO ji = 1, jpi
110                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
111               ENDDO
112            END DO
113            CALL iom_put( 'zomsf', z4d1 * rc_sv )
114         ENDIF
115         IF(  iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.   &
116            & iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
117            ! define fields multiplied by scalar
118            zmask(:,:,:) = 0._wp
119            zts(:,:,:,:) = 0._wp
120            DO_3D_10_11( 1, jpkm1 )
121               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
122               zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
123               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
124               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
125            END_3D
126         ENDIF
127         IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
128            DO jn = 1, nptr
129               sjk(:,:,jn) = ptr_sjk( zmask(:,:,:), btmsk(:,:,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) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
134               zs_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
135               v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk34(:,:,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) = ptr_sj( zmask(:,:,:), btmsk(:,:,jn) )
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(:) = ptr_sj( pvtr(:,:,:), btmsk34(:,:,jn) )
164               ztsum(:) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
165               zssum(:) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
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      ELSE
187         !
188         zmask(:,:,:) = 0._wp
189         zts(:,:,:,:) = 0._wp
190         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface
191            DO_3D_11_11( 1, jpkm1 )
192               zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm)
193               zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
194               zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc
195               zts(ji,jj,jk,jp_sal) = ts(ji,jj,jk,jp_sal,Kmm) * zsfc
196            END_3D
197            !
198            DO jn = 1, nptr
199               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
200               z4d1(:,:,:,jn) = zmask(:,:,:)
201            ENDDO
202            CALL iom_put( 'zosrf', z4d1 )
203            !
204            DO jn = 1, nptr
205               z4d2(1,:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) &
206                  &            / MAX( z4d1(1,:,:,jn), 10.e-15 )
207               DO ji = 1, jpi
208                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
209               ENDDO
210            ENDDO
211            CALL iom_put( 'zotem', z4d2 )
212            !
213            DO jn = 1, nptr
214               z4d2(1,:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) &
215                  &            / MAX( z4d1(1,:,:,jn), 10.e-15 )
216               DO ji = 1, jpi
217                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
218               ENDDO
219            ENDDO
220            CALL iom_put( 'zosal', z4d2 )
221            !
222         ENDIF
223         !
224         !                                ! Advective and diffusive heat and salt transport
225         IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN 
226            !
227            DO jn = 1, nptr
228               z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
229               DO ji = 1, jpi
230                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
231               ENDDO
232            ENDDO
233            CALL iom_put( 'sophtadv', z3dtr )
234            DO jn = 1, nptr
235               z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
236               DO ji = 1, jpi
237                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
238               ENDDO
239            ENDDO
240            CALL iom_put( 'sopstadv', z3dtr )
241         ENDIF
242         !
243         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN 
244            !
245            DO jn = 1, nptr
246               z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
247               DO ji = 1, jpi
248                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
249               ENDDO
250            ENDDO
251            CALL iom_put( 'sophtldf', z3dtr )
252            DO jn = 1, nptr
253               z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
254               DO ji = 1, jpi
255                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
256               ENDDO
257            ENDDO
258            CALL iom_put( 'sopstldf', z3dtr )
259         ENDIF
260         !
261         IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN 
262            !
263            DO jn = 1, nptr
264               z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
265               DO ji = 1, jpi
266                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
267               ENDDO
268            ENDDO
269            CALL iom_put( 'sophteiv', z3dtr )
270            DO jn = 1, nptr
271               z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
272               DO ji = 1, jpi
273                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
274               ENDDO
275            ENDDO
276            CALL iom_put( 'sopsteiv', z3dtr )
277         ENDIF
278         !
279         IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
280            zts(:,:,:,:) = 0._wp
281            DO_3D_10_11( 1, jpkm1 )
282               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
283               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
284               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
285            END_3D
286             CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) )
287             CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) )
288             DO jn = 1, nptr
289                z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
290                DO ji = 1, jpi
291                   z3dtr(ji,:,jn) = z3dtr(1,:,jn)
292                ENDDO
293             ENDDO
294             CALL iom_put( 'sophtvtr', z3dtr )
295             DO jn = 1, nptr
296               z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
297               DO ji = 1, jpi
298                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
299               ENDDO
300            ENDDO
301            CALL iom_put( 'sopstvtr', z3dtr )
302         ENDIF
303         !
304         IF( iom_use( 'uocetr_vsum_cumul' ) ) THEN
305            CALL iom_get_var(  'uocetr_vsum_op', z2d ) ! get uocetr_vsum_op from xml
306            z2d(:,:) = ptr_ci_2d( z2d(:,:) ) 
307            CALL iom_put( 'uocetr_vsum_cumul', z2d )
308         ENDIF
309         !
310      ENDIF
311      !
312      IF( ln_timing )   CALL timing_stop('dia_ptr')
313      !
314   END SUBROUTINE dia_ptr
315
316
317   SUBROUTINE dia_ptr_init
318      !!----------------------------------------------------------------------
319      !!                  ***  ROUTINE dia_ptr_init  ***
320      !!                   
321      !! ** Purpose :   Initialization, namelist read
322      !!----------------------------------------------------------------------
323      INTEGER ::  inum, jn           ! local integers
324      !!
325      REAL(wp), DIMENSION(jpi,jpj) :: zmsk
326      !!----------------------------------------------------------------------
327
328      l_diaptr = .FALSE.
329      IF(   iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  &
330         &  iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  &
331         &  iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  &
332         &  iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  & 
333         &  iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  &
334         &  iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) )  l_diaptr  = .TRUE.
335
336 
337      IF(lwp) THEN                     ! Control print
338         WRITE(numout,*)
339         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
340         WRITE(numout,*) '~~~~~~~~~~~~'
341         WRITE(numout,*) '   Namelist namptr : set ptr parameters'
342         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr
343      ENDIF
344
345      IF( l_diaptr ) THEN 
346         !
347         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
348
349         rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt
350         rc_ggram = rc_ggram * rau0              ! conversion from m3/s to Gg/s
351
352         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
353
354         btmsk(:,:,1) = tmask_i(:,:)                 
355         CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  )
356         CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
357         CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
358         CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
359         CALL iom_close( inum )
360         btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
361         DO jn = 2, nptr
362            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
363         END DO
364         ! JD : modification so that overturning streamfunction is available in Atlantic at 34S to compare with observations
365         WHERE( gphit(:,:)*tmask_i(:,:) < -34._wp)
366           zmsk(:,:) = 0._wp      ! mask out Southern Ocean
367         ELSE WHERE                 
368           zmsk(:,:) = ssmask(:,:)
369         END WHERE
370         btmsk34(:,:,1) = btmsk(:,:,1)                 
371         DO jn = 2, nptr
372            btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:)               ! interior domain only
373         ENDDO
374
375         ! Initialise arrays to zero because diatpr is called before they are first calculated
376         ! Note that this means diagnostics will not be exactly correct when model run is restarted.
377         hstr_adv(:,:,:) = 0._wp           
378         hstr_ldf(:,:,:) = 0._wp           
379         hstr_eiv(:,:,:) = 0._wp           
380         hstr_ove(:,:,:) = 0._wp           
381         hstr_btr(:,:,:) = 0._wp           !
382         hstr_vtr(:,:,:) = 0._wp           !
383         !
384         ll_init = .FALSE.
385         !
386      ENDIF 
387      !
388   END SUBROUTINE dia_ptr_init
389
390
391   SUBROUTINE dia_ptr_hst( ktra, cptr, pvflx ) 
392      !!----------------------------------------------------------------------
393      !!                    ***  ROUTINE dia_ptr_hst ***
394      !!----------------------------------------------------------------------
395      !! Wrapper for heat and salt transport calculations to calculate them for each basin
396      !! Called from all advection and/or diffusion routines
397      !!----------------------------------------------------------------------
398      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
399      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
400      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pvflx   ! 3D input array of advection/diffusion
401      INTEGER                                        :: jn    !
402
403      !
404      IF( cptr == 'adv' ) THEN
405         IF( ktra == jp_tem )  THEN
406             DO jn = 1, nptr
407                hstr_adv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
408             ENDDO
409         ENDIF
410         IF( ktra == jp_sal )  THEN
411             DO jn = 1, nptr
412                hstr_adv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
413             ENDDO
414         ENDIF
415      ENDIF
416      !
417      IF( cptr == 'ldf' ) THEN
418         IF( ktra == jp_tem )  THEN
419             DO jn = 1, nptr
420                hstr_ldf(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
421             ENDDO
422         ENDIF
423         IF( ktra == jp_sal )  THEN
424             DO jn = 1, nptr
425                hstr_ldf(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
426             ENDDO
427         ENDIF
428      ENDIF
429      !
430      IF( cptr == 'eiv' ) THEN
431         IF( ktra == jp_tem )  THEN
432             DO jn = 1, nptr
433                hstr_eiv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
434             ENDDO
435         ENDIF
436         IF( ktra == jp_sal )  THEN
437             DO jn = 1, nptr
438                hstr_eiv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
439             ENDDO
440         ENDIF
441      ENDIF
442      !
443      IF( cptr == 'vtr' ) THEN
444         IF( ktra == jp_tem )  THEN
445             DO jn = 1, nptr
446                hstr_vtr(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
447             ENDDO
448         ENDIF
449         IF( ktra == jp_sal )  THEN
450             DO jn = 1, nptr
451                hstr_vtr(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
452             ENDDO
453         ENDIF
454      ENDIF
455      !
456   END SUBROUTINE dia_ptr_hst
457
458
459   FUNCTION dia_ptr_alloc()
460      !!----------------------------------------------------------------------
461      !!                    ***  ROUTINE dia_ptr_alloc  ***
462      !!----------------------------------------------------------------------
463      INTEGER               ::   dia_ptr_alloc   ! return value
464      INTEGER, DIMENSION(3) ::   ierr
465      !!----------------------------------------------------------------------
466      ierr(:) = 0
467      !
468      IF( .NOT. ALLOCATED( btmsk ) ) THEN
469         ALLOCATE( btmsk(jpi,jpj,nptr)    , btmsk34(jpi,jpj,nptr),   &
470            &      hstr_adv(jpj,jpts,nptr), hstr_eiv(jpj,jpts,nptr), &
471            &      hstr_ove(jpj,jpts,nptr), hstr_btr(jpj,jpts,nptr), &
472            &      hstr_ldf(jpj,jpts,nptr), hstr_vtr(jpj,jpts,nptr), STAT=ierr(1)  )
473            !
474         ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2))
475         !
476         dia_ptr_alloc = MAXVAL( ierr )
477         CALL mpp_sum( 'diaptr', dia_ptr_alloc )
478      ENDIF
479      !
480   END FUNCTION dia_ptr_alloc
481
482
483   FUNCTION ptr_sj_3d( pvflx, pmsk )   RESULT ( p_fval )
484      !!----------------------------------------------------------------------
485      !!                    ***  ROUTINE ptr_sj_3d  ***
486      !!
487      !! ** Purpose :   i-k sum computation of a j-flux array
488      !!
489      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
490      !!              pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
491      !!
492      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
493      !!----------------------------------------------------------------------
494      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)  ::   pvflx  ! mask flux array at V-point
495      REAL(wp), INTENT(in), DIMENSION(jpi,jpj)      ::   pmsk   ! Optional 2D basin mask
496      !
497      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
498      INTEGER                  ::   ijpj         ! ???
499      REAL(wp), POINTER, DIMENSION(:) :: p_fval  ! function value
500      !!--------------------------------------------------------------------
501      !
502      p_fval => p_fval1d
503
504      ijpj = jpj
505      p_fval(:) = 0._wp
506      DO_3D_00_00( 1, jpkm1 )
507         p_fval(jj) = p_fval(jj) + pvflx(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
508      END_3D
509#if defined key_mpp_mpi
510      CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl)
511#endif
512      !
513   END FUNCTION ptr_sj_3d
514
515
516   FUNCTION ptr_sj_2d( pvflx, pmsk )   RESULT ( p_fval )
517      !!----------------------------------------------------------------------
518      !!                    ***  ROUTINE ptr_sj_2d  ***
519      !!
520      !! ** Purpose :   "zonal" and vertical sum computation of a j-flux array
521      !!
522      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
523      !!      pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
524      !!
525      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
526      !!----------------------------------------------------------------------
527      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pvflx  ! mask flux array at V-point
528      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
529      !
530      INTEGER                  ::   ji,jj       ! dummy loop arguments
531      INTEGER                  ::   ijpj        ! ???
532      REAL(wp), POINTER, DIMENSION(:) :: p_fval ! function value
533      !!--------------------------------------------------------------------
534      !
535      p_fval => p_fval1d
536
537      ijpj = jpj
538      p_fval(:) = 0._wp
539      DO_2D_00_00
540         p_fval(jj) = p_fval(jj) + pvflx(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj)
541      END_2D
542#if defined key_mpp_mpi
543      CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl )
544#endif
545      !
546   END FUNCTION ptr_sj_2d
547
548   FUNCTION ptr_ci_2d( pva )   RESULT ( p_fval )
549      !!----------------------------------------------------------------------
550      !!                    ***  ROUTINE ptr_ci_2d  ***
551      !!
552      !! ** Purpose :   "meridional" cumulated sum computation of a j-flux array
553      !!
554      !! ** Method  : - j cumulated sum of pva using the interior 2D vmask (umask_i).
555      !!
556      !! ** Action  : - p_fval: j-cumulated sum of pva
557      !!----------------------------------------------------------------------
558      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)  ::   pva   ! mask flux array at V-point
559      !
560      INTEGER                  ::   ji,jj,jc       ! dummy loop arguments
561      INTEGER                  ::   ijpj        ! ???
562      REAL(wp), DIMENSION(jpi,jpj) :: p_fval ! function value
563      !!--------------------------------------------------------------------
564      !
565      ijpj = jpj  ! ???
566      p_fval(:,:) = 0._wp
567      DO jc = 1, jpnj ! looping over all processors in j axis
568         DO_2D_00_00
569            p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj)
570         END_2D
571         CALL lbc_lnk( 'diaptr', p_fval, 'U', -1. )
572      END DO
573      !
574   END FUNCTION ptr_ci_2d
575
576
577
578   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
579      !!----------------------------------------------------------------------
580      !!                    ***  ROUTINE ptr_sjk  ***
581      !!
582      !! ** Purpose :   i-sum computation of an array
583      !!
584      !! ** Method  : - i-sum of field using the interior 2D vmask (pmsk).
585      !!
586      !! ** Action  : - p_fval: i-sum of masked field
587      !!----------------------------------------------------------------------
588      !!
589      IMPLICIT none
590      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! mask flux array at V-point
591      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)     ::   pmsk   ! Optional 2D basin mask
592      !!
593      INTEGER                           :: ji, jj, jk ! dummy loop arguments
594      REAL(wp), POINTER, DIMENSION(:,:) :: p_fval     ! return function value
595#if defined key_mpp_mpi
596      INTEGER, DIMENSION(1) ::   ish
597      INTEGER, DIMENSION(2) ::   ish2
598      INTEGER               ::   ijpjjpk
599      REAL(wp), DIMENSION(jpj*jpk) ::   zwork    ! mask flux array at V-point
600#endif
601      !!--------------------------------------------------------------------
602      !
603      p_fval => p_fval2d
604
605      p_fval(:,:) = 0._wp
606      !
607      DO_3D_00_00( 1, jpkm1 )
608         p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
609      END_3D
610      !
611#if defined key_mpp_mpi
612      ijpjjpk = jpj*jpk
613      ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk
614      zwork(1:ijpjjpk) = RESHAPE( p_fval, ish )
615      CALL mpp_sum( 'diaptr', zwork, ijpjjpk, ncomm_znl )
616      p_fval(:,:) = RESHAPE( zwork, ish2 )
617#endif
618      !
619   END FUNCTION ptr_sjk
620
621
622   !!======================================================================
623END MODULE diaptr
Note: See TracBrowser for help on using the repository browser.