source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

File size: 36.3 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   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   dia_ptr      : Poleward Transport Diagnostics module
16   !!   dia_ptr_init : Initialization, namelist read
17   !!   ptr_sjk      : "zonal" mean computation of a field - tracer or flux array
18   !!   ptr_sj       : "zonal" and vertical sum computation of a "meridional" flux array
19   !!                   (Generic interface to ptr_sj_3d, ptr_sj_2d)
20   !!----------------------------------------------------------------------
21   USE oce              ! ocean dynamics and active tracers
22   USE dom_oce          ! ocean space and time domain
23   USE phycst           ! physical constants
24   USE ldftra_oce 
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   USE yomhook, ONLY: lhook, dr_hook
32   USE parkind1, ONLY: jprb, jpim
33
34   IMPLICIT NONE
35   PRIVATE
36
37   INTERFACE ptr_sj
38      MODULE PROCEDURE ptr_sj_3d, ptr_sj_2d
39   END INTERFACE
40
41   PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines
42   PUBLIC   ptr_sjk        !
43   PUBLIC   dia_ptr_init   ! call in step module
44   PUBLIC   dia_ptr        ! call in step module
45   PUBLIC   dia_ptr_ohst_components        ! called from tra_ldf/tra_adv routines
46
47   !                                  !!** namelist  namptr  **
48   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_adv, htr_ldf, htr_eiv, htr_vt   !: Heat TRansports (adv, diff, Bolus.)
49   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   str_adv, str_ldf, str_eiv, str_vs   !: Salt TRansports (adv, diff, Bolus.)
50   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_ove, str_ove   !: heat Salt TRansports ( overturn.)
51   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_btr, str_btr   !: heat Salt TRansports ( barotropic )
52
53   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F)
54   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation
55   INTEGER, PUBLIC ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T)
56
57   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
58   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rau0 x Cp)
59   REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg
60
61   CHARACTER(len=3), ALLOCATABLE, SAVE, DIMENSION(:)     :: clsubb
62   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
63   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   :: btm30   ! mask out Southern Ocean (=0 south of 30°S)
64
65   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)     :: p_fval1d
66   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: p_fval2d
67
68
69   !! * Substitutions
70#  include "domzgr_substitute.h90"
71#  include "vectopt_loop_substitute.h90"
72   !!----------------------------------------------------------------------
73   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
74   !! $Id$
75   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
76   !!----------------------------------------------------------------------
77CONTAINS
78
79   SUBROUTINE dia_ptr( pvtr )
80      !!----------------------------------------------------------------------
81      !!                  ***  ROUTINE dia_ptr  ***
82      !!----------------------------------------------------------------------
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      REAL(wp) ::   zsfc,zvfc               ! local scalar
87      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
88      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace
89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace
90      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace
91      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace
92      REAL(wp), DIMENSION(jpj)     ::  vsum   ! 1D workspace
93      REAL(wp), DIMENSION(jpj,jpts)     ::  tssum   ! 1D workspace
94 
95      !
96      !overturning calculation
97      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   sjk  , r1_sjk ! i-mean i-k-surface and its inverse
98      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   v_msf, sn_jk  , tn_jk ! i-mean T and S, j-Stream-Function
99
100
101      CHARACTER( len = 12 )  :: cl1
102      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
103      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
104      REAL(KIND=jprb)               :: zhook_handle
105
106      CHARACTER(LEN=*), PARAMETER :: RoutineName='DIA_PTR'
107
108      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
109
110      !!----------------------------------------------------------------------
111      !
112      IF( nn_timing == 1 )   CALL timing_start('dia_ptr')
113
114      !
115      z2d(:,:) = 0._wp
116      z3d(:,:,:) = 0._wp
117      IF( PRESENT( pvtr ) ) THEN
118         IF( iom_use("zomsfglo") ) THEN    ! effective MSF
119            z3d(1,:,:) = ptr_sjk( pvtr(:,:,:) )  ! zonal cumulative effective transport
120            DO jk = jpkm1,1,-1                   !Integrate from bottom up to get
121              z3d(1,:,jk) = z3d(1,:,jk+1) - z3d(1,:,jk)   ! effective j-Stream-Function (MSF)
122            END DO
123            DO ji = 1, jpi
124               z3d(ji,:,:) = z3d(1,:,:)
125            ENDDO
126            cl1 = TRIM('zomsf'//clsubb(1) )
127            CALL iom_put( cl1, z3d * rc_sv )
128            DO jn = 2, nptr                                    ! by sub-basins
129               z3d(1,:,:) =  ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn) ) 
130               DO jk = jpkm1,1,-1
131                  z3d(1,:,jk) = z3d(1,:,jk+1) - z3d(1,:,jk)    ! effective j-Stream-Function (MSF)
132               END DO
133               DO ji = 1, jpi
134                  z3d(ji,:,:) = z3d(1,:,:)
135               ENDDO
136               cl1 = TRIM('zomsf'//clsubb(jn) )
137               CALL iom_put( cl1, z3d * rc_sv )
138            END DO
139         ENDIF
140         IF( iom_use("sopstove") .OR. iom_use("sophtove") .OR. iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN
141            ! define fields multiplied by scalar
142            zmask(:,:,:) = 0._wp
143            zts(:,:,:,:) = 0._wp
144            DO jk = 1, jpkm1
145               DO jj = 1, jpjm1
146                  DO ji = 1, jpi
147                     zvfc = e1v(ji,jj) * fse3v(ji,jj,jk)
148                     zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
149                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
150                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
151                  ENDDO
152               ENDDO
153             ENDDO
154         ENDIF
155         IF( iom_use("sopstove") .OR. iom_use("sophtove") ) THEN
156             sjk(:,:,1) = ptr_sjk( zmask(:,:,:), btmsk(:,:,1) )
157             r1_sjk(:,:,1) = 0._wp
158             WHERE( sjk(:,:,1) /= 0._wp )   r1_sjk(:,:,1) = 1._wp / sjk(:,:,1)
159
160             ! i-mean T and S, j-Stream-Function, global
161             tn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_tem) ) * r1_sjk(:,:,1)
162             sn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_sal) ) * r1_sjk(:,:,1)
163             v_msf(:,:,1) = ptr_sjk( pvtr(:,:,:) )
164
165             htr_ove(:,1) = SUM( v_msf(:,:,1)*tn_jk(:,:,1) ,2 )
166             str_ove(:,1) = SUM( v_msf(:,:,1)*sn_jk(:,:,1) ,2 )
167
168             z2d(1,:) = htr_ove(:,1) * rc_pwatt        !  (conversion in PW)
169             DO ji = 1, jpi
170               z2d(ji,:) = z2d(1,:)
171             ENDDO
172             cl1 = 'sophtove'
173             CALL iom_put( TRIM(cl1), z2d )
174             z2d(1,:) = str_ove(:,1) * rc_ggram        !  (conversion in Gg)
175             DO ji = 1, jpi
176               z2d(ji,:) = z2d(1,:)
177             ENDDO
178             cl1 = 'sopstove'
179             CALL iom_put( TRIM(cl1), z2d )
180             IF( ln_subbas ) THEN
181                DO jn = 2, nptr
182                    sjk(:,:,jn) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
183                    r1_sjk(:,:,jn) = 0._wp
184                    WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
185
186                    ! i-mean T and S, j-Stream-Function, basin
187                    tn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
188                    sn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
189                    v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn) ) 
190                    htr_ove(:,jn) = SUM( v_msf(:,:,jn)*tn_jk(:,:,jn) ,2 )
191                    str_ove(:,jn) = SUM( v_msf(:,:,jn)*sn_jk(:,:,jn) ,2 )
192
193                    z2d(1,:) = htr_ove(:,jn) * rc_pwatt !  (conversion in PW)
194                    DO ji = 1, jpi
195                        z2d(ji,:) = z2d(1,:)
196                    ENDDO
197                    cl1 = TRIM('sophtove_'//clsubb(jn))
198                    CALL iom_put( cl1, z2d )
199                    z2d(1,:) = str_ove(:,jn) * rc_ggram        ! (conversion in Gg)
200                    DO ji = 1, jpi
201                        z2d(ji,:) = z2d(1,:)
202                    ENDDO
203                    cl1 = TRIM('sopstove_'//clsubb(jn))
204                    CALL iom_put( cl1, z2d )
205                END DO
206             ENDIF
207         ENDIF
208         IF( iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN
209         ! Calculate barotropic heat and salt transport here
210             sjk(:,1,1) = ptr_sj( zmask(:,:,:), btmsk(:,:,1) )
211             r1_sjk(:,1,1) = 0._wp
212             WHERE( sjk(:,1,1) /= 0._wp )   r1_sjk(:,1,1) = 1._wp / sjk(:,1,1)
213           
214            vsum = ptr_sj( pvtr(:,:,:), btmsk(:,:,1))
215            tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,1) )
216            tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,1) )
217            htr_btr(:,1) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,1)
218            str_btr(:,1) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,1)
219            z2d(1,:) = htr_btr(:,1) * rc_pwatt        !  (conversion in PW)
220            DO ji = 2, jpi
221               z2d(ji,:) = z2d(1,:)
222            ENDDO
223            cl1 = 'sophtbtr'
224            CALL iom_put( TRIM(cl1), z2d )
225            z2d(1,:) = str_btr(:,1) * rc_ggram        !  (conversion in Gg)
226            DO ji = 2, jpi
227              z2d(ji,:) = z2d(1,:)
228            ENDDO
229            cl1 = 'sopstbtr'
230            CALL iom_put( TRIM(cl1), z2d )
231            IF( ln_subbas ) THEN
232                DO jn = 2, nptr
233                    sjk(:,1,jn) = ptr_sj( zmask(:,:,:), btmsk(:,:,jn) )
234                    r1_sjk(:,1,jn) = 0._wp
235                    WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
236                    vsum = ptr_sj( pvtr(:,:,:), btmsk(:,:,jn))
237                    tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
238                    tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
239                    htr_btr(:,jn) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,jn)
240                    str_btr(:,jn) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,jn)
241                    z2d(1,:) = htr_btr(:,jn) * rc_pwatt !  (conversion in PW)
242                    DO ji = 1, jpi
243                        z2d(ji,:) = z2d(1,:)
244                    ENDDO
245                    cl1 = TRIM('sophtbtr_'//clsubb(jn))
246                    CALL iom_put( cl1, z2d )
247                    z2d(1,:) = str_btr(:,jn) * rc_ggram        ! (conversion in Gg)
248                    DO ji = 1, jpi
249                        z2d(ji,:) = z2d(1,:)
250                    ENDDO
251                    cl1 = TRIM('sopstbtr_'//clsubb(jn))
252                    CALL iom_put( cl1, z2d )
253               ENDDO
254            ENDIF !ln_subbas
255         ENDIF !iom_use("sopstbtr....)
256         !
257      ELSE
258         !
259         IF( iom_use("zotemglo") ) THEN    ! i-mean i-k-surface
260            zmask(:,:,:) = 0._wp
261            zts(:,:,:,:) = 0._wp
262            DO jk = 1, jpkm1
263               DO jj = 1, jpj
264                  DO ji = 1, jpi
265                     zsfc = e1t(ji,jj) * fse3t(ji,jj,jk)
266                     zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
267                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc
268                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc
269                  ENDDO
270               ENDDO
271            ENDDO
272            DO jn = 1, nptr
273               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
274               cl1 = TRIM('zosrf'//clsubb(jn) )
275               CALL iom_put( cl1, zmask )
276               !
277               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) &
278                  &            / MAX( zmask(1,:,:), 10.e-15 )
279               DO ji = 1, jpi
280                  z3d(ji,:,:) = z3d(1,:,:)
281               ENDDO
282               cl1 = TRIM('zotem'//clsubb(jn) )
283               CALL iom_put( cl1, z3d )
284               !
285               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) &
286                  &            / MAX( zmask(1,:,:), 10.e-15 )
287               DO ji = 1, jpi
288                  z3d(ji,:,:) = z3d(1,:,:)
289               ENDDO
290               cl1 = TRIM('zosal'//clsubb(jn) )
291               CALL iom_put( cl1, z3d )
292            END DO
293         ENDIF
294         !
295         !                                ! Advective and diffusive heat and salt transport
296         IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN   
297            z2d(1,:) = htr_adv(:,1) * rc_pwatt        !  (conversion in PW)
298            DO ji = 1, jpi
299               z2d(ji,:) = z2d(1,:)
300            ENDDO
301            cl1 = 'sophtadv'                 
302            CALL iom_put( TRIM(cl1), z2d )
303            z2d(1,:) = str_adv(:,1) * rc_ggram        ! (conversion in Gg)
304            DO ji = 1, jpi
305               z2d(ji,:) = z2d(1,:)
306            ENDDO
307            cl1 = 'sopstadv'
308            CALL iom_put( TRIM(cl1), z2d )
309            IF( ln_subbas ) THEN
310              DO jn=2,nptr
311               z2d(1,:) = htr_adv(:,jn) * rc_pwatt        !  (conversion in PW)
312               DO ji = 1, jpi
313                 z2d(ji,:) = z2d(1,:)
314               ENDDO
315               cl1 = TRIM('sophtadv_'//clsubb(jn))                 
316               CALL iom_put( cl1, z2d )
317               z2d(1,:) = str_adv(:,jn) * rc_ggram        ! (conversion in Gg)
318               DO ji = 1, jpi
319                  z2d(ji,:) = z2d(1,:)
320               ENDDO
321               cl1 = TRIM('sopstadv_'//clsubb(jn))                 
322               CALL iom_put( cl1, z2d )             
323              ENDDO
324            ENDIF
325         ENDIF
326         !
327         IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN   
328            z2d(1,:) = htr_ldf(:,1) * rc_pwatt        !  (conversion in PW)
329            DO ji = 1, jpi
330               z2d(ji,:) = z2d(1,:)
331            ENDDO
332            cl1 = 'sophtldf'
333            CALL iom_put( TRIM(cl1), z2d )
334            z2d(1,:) = str_ldf(:,1) * rc_ggram        !  (conversion in Gg)
335            DO ji = 1, jpi
336               z2d(ji,:) = z2d(1,:)
337            ENDDO
338            cl1 = 'sopstldf'
339            CALL iom_put( TRIM(cl1), z2d )
340            IF( ln_subbas ) THEN
341              DO jn=2,nptr
342               z2d(1,:) = htr_ldf(:,jn) * rc_pwatt        !  (conversion in PW)
343               DO ji = 1, jpi
344                 z2d(ji,:) = z2d(1,:)
345               ENDDO
346               cl1 = TRIM('sophtldf_'//clsubb(jn))                 
347               CALL iom_put( cl1, z2d )
348               z2d(1,:) = str_ldf(:,jn) * rc_ggram        ! (conversion in Gg)
349               DO ji = 1, jpi
350                  z2d(ji,:) = z2d(1,:)
351               ENDDO
352               cl1 = TRIM('sopstldf_'//clsubb(jn))                 
353               CALL iom_put( cl1, z2d )             
354              ENDDO
355            ENDIF
356         ENDIF
357
358         IF( iom_use("sopht_vt") .OR. iom_use("sopst_vs") ) THEN   
359            z2d(1,:) = htr_vt(:,1) * rc_pwatt        !  (conversion in PW)
360            DO ji = 1, jpi
361               z2d(ji,:) = z2d(1,:)
362            ENDDO
363            cl1 = 'sopht_vt'
364            CALL iom_put( TRIM(cl1), z2d )
365            z2d(1,:) = str_vs(:,1) * rc_ggram        !  (conversion in Gg)
366            DO ji = 1, jpi
367               z2d(ji,:) = z2d(1,:)
368            ENDDO
369            cl1 = 'sopst_vs'
370            CALL iom_put( TRIM(cl1), z2d )
371            IF( ln_subbas ) THEN
372              DO jn=2,nptr
373               z2d(1,:) = htr_vt(:,jn) * rc_pwatt        !  (conversion in PW)
374               DO ji = 1, jpi
375                 z2d(ji,:) = z2d(1,:)
376               ENDDO
377               cl1 = TRIM('sopht_vt_'//clsubb(jn))                 
378               CALL iom_put( cl1, z2d )
379               z2d(1,:) = str_vs(:,jn) * rc_ggram        ! (conversion in Gg)
380               DO ji = 1, jpi
381                  z2d(ji,:) = z2d(1,:)
382               ENDDO
383               cl1 = TRIM('sopst_vs_'//clsubb(jn))                 
384               CALL iom_put( cl1, z2d )             
385              ENDDO
386            ENDIF
387         ENDIF
388
389#ifdef key_diaeiv
390         IF(lk_traldf_eiv) THEN
391            IF( iom_use("sophteiv") .OR. iom_use("sopsteiv") ) THEN
392               z2d(1,:) = htr_eiv(:,1) * rc_pwatt        !  (conversion in PW)
393               DO ji = 1, jpi
394                  z2d(ji,:) = z2d(1,:)
395               ENDDO
396               cl1 = 'sophteiv'
397               CALL iom_put( TRIM(cl1), z2d )
398               z2d(1,:) = str_eiv(:,1) * rc_ggram        !  (conversion in Gg)
399               DO ji = 1, jpi
400                  z2d(ji,:) = z2d(1,:)
401               ENDDO
402               cl1 = 'sopsteiv'
403               CALL iom_put( TRIM(cl1), z2d )
404               IF( ln_subbas ) THEN
405                  DO jn=2,nptr
406                     z2d(1,:) = htr_eiv(:,jn) * rc_pwatt        !  (conversion in PW)
407                     DO ji = 1, jpi
408                        z2d(ji,:) = z2d(1,:)
409                     ENDDO
410                     cl1 = TRIM('sophteiv_'//clsubb(jn))                 
411                     CALL iom_put( cl1, z2d )
412                     z2d(1,:) = str_eiv(:,jn) * rc_ggram        ! (conversion in Gg)
413                     DO ji = 1, jpi
414                        z2d(ji,:) = z2d(1,:)
415                     ENDDO
416                     cl1 = TRIM('sopsteiv_'//clsubb(jn)) 
417                     CALL iom_put( cl1, z2d )             
418                  ENDDO
419               ENDIF
420            ENDIF
421            IF( iom_use("zomsfeivglo") ) THEN
422               DO jk=1,jpk
423                  DO jj=1,jpj
424                     DO ji=1,jpi
425                        zvn(ji,jj,jk) = v_eiv(ji,jj,jk) * fse3v(ji,jj,jk) * e1v(ji,jj)
426                     ENDDO
427                  ENDDO
428               ENDDO
429               z3d(1,:,:) = ptr_sjk( zvn(:,:,:) )  ! zonal cumulative effective transport
430               DO jk = jpkm1,1,-1
431                 z3d(1,:,jk) = z3d(1,:,jk+1) - z3d(1,:,jk)   ! effective j-Stream-Function (MSF)
432               END DO
433               DO ji = 1, jpi
434                  z3d(ji,:,:) = z3d(1,:,:)
435               ENDDO
436               cl1 = TRIM('zomsfeiv'//clsubb(1) )
437               CALL iom_put( cl1, z3d * rc_sv )
438               IF( ln_subbas ) THEN
439                  DO jn = 2, nptr                                    ! by sub-basins
440                     z3d(1,:,:) =  ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) ) 
441                     DO jk = jpkm1,1,-1
442                        z3d(1,:,jk) = z3d(1,:,jk+1) - z3d(1,:,jk)    ! effective j-Stream-Function (MSF)
443                     END DO
444                     DO ji = 1, jpi
445                        z3d(ji,:,:) = z3d(1,:,:)
446                     ENDDO
447                     cl1 = TRIM('zomsfeiv'//clsubb(jn) )
448                     CALL iom_put( cl1, z3d * rc_sv )
449                  END DO
450               ENDIF
451            ENDIF
452         ENDIF
453#endif
454         !
455      ENDIF
456      !
457      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr')
458      !
459      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
460   END SUBROUTINE dia_ptr
461
462
463   SUBROUTINE dia_ptr_init
464      !!----------------------------------------------------------------------
465      !!                  ***  ROUTINE dia_ptr_init  ***
466      !!                   
467      !! ** Purpose :   Initialization, namelist read
468      !!----------------------------------------------------------------------
469      INTEGER ::  jn           ! local integers
470      INTEGER ::  inum, ierr   ! local integers
471      INTEGER ::  ios          ! Local integer output status for namelist read
472      !!
473      NAMELIST/namptr/ ln_diaptr, ln_subbas
474      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
475      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
476      REAL(KIND=jprb)               :: zhook_handle
477
478      CHARACTER(LEN=*), PARAMETER :: RoutineName='DIA_PTR_INIT'
479
480      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
481
482      !!----------------------------------------------------------------------
483
484      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
485      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901)
486901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist', lwp )
487
488      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport
489      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 )
490902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp )
491      IF(lwm) WRITE ( numond, namptr )
492
493      IF(lwp) THEN                     ! Control print
494         WRITE(numout,*)
495         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
496         WRITE(numout,*) '~~~~~~~~~~~~'
497         WRITE(numout,*) '   Namelist namptr : set ptr parameters'
498         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr
499         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas
500      ENDIF
501
502      IF( ln_diaptr ) THEN 
503         !
504         IF( ln_subbas ) THEN
505            nptr = 5            ! Global, Atlantic, Pacific, Indian, Indo-Pacific
506            ALLOCATE( clsubb(nptr) )
507            clsubb(1) = 'glo' ;  clsubb(2) = 'atl'  ;  clsubb(3) = 'pac'  ;  clsubb(4) = 'ind'  ;  clsubb(5) = 'ipc'
508         ELSE               
509            nptr = 1       ! Global only
510            ALLOCATE( clsubb(nptr) )
511            clsubb(1) = 'glo' 
512         ENDIF
513
514         !                                      ! allocate dia_ptr arrays
515         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
516
517         rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt
518
519         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
520
521         IF( ln_subbas ) THEN                ! load sub-basin mask
522            CALL iom_open( 'subbasins', inum,  ldstop = .TRUE.  )
523            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
524            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
525            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
526            CALL iom_close( inum )
527            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
528            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
529            ELSE WHERE                     ;   btm30(:,:) = ssmask(:,:)
530            END WHERE
531         ENDIF
532   
533         btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
534     
535         DO jn = 1, nptr
536            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
537         END DO
538
539         ! Initialise arrays to zero because diatpr is called before they are first calculated
540         ! Note that this means diagnostics will not be exactly correct when model run is restarted.
541         htr_adv(:,:) = 0._wp  ;  str_adv(:,:) =  0._wp 
542         htr_ldf(:,:) = 0._wp  ;  str_ldf(:,:) =  0._wp 
543         htr_eiv(:,:) = 0._wp  ;  str_eiv(:,:) =  0._wp 
544         htr_vt(:,:) = 0._wp  ;   str_vs(:,:) =  0._wp
545         htr_ove(:,:) = 0._wp  ;   str_ove(:,:) =  0._wp
546         htr_btr(:,:) = 0._wp  ;   str_btr(:,:) =  0._wp
547         !
548      ENDIF 
549      !
550      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
551   END SUBROUTINE dia_ptr_init
552
553   SUBROUTINE dia_ptr_ohst_components( ktra, cptr, pva ) 
554      !!----------------------------------------------------------------------
555      !!                    ***  ROUTINE dia_ptr_ohst_components  ***
556      !!----------------------------------------------------------------------
557      !! Wrapper for heat and salt transport calculations to calculate them for each basin
558      !! Called from all advection and/or diffusion routines
559      !!----------------------------------------------------------------------
560      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
561      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
562      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
563      INTEGER                                        :: jn    !
564      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
565      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
566      REAL(KIND=jprb)               :: zhook_handle
567
568      CHARACTER(LEN=*), PARAMETER :: RoutineName='DIA_PTR_OHST_COMPONENTS'
569
570      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
571
572
573      IF( cptr == 'adv' ) THEN
574         IF( ktra == jp_tem )  htr_adv(:,1) = ptr_sj( pva(:,:,:) )
575         IF( ktra == jp_sal )  str_adv(:,1) = ptr_sj( pva(:,:,:) )
576      ENDIF
577      IF( cptr == 'ldf' ) THEN
578         IF( ktra == jp_tem )  htr_ldf(:,1) = ptr_sj( pva(:,:,:) )
579         IF( ktra == jp_sal )  str_ldf(:,1) = ptr_sj( pva(:,:,:) )
580      ENDIF
581      IF( cptr == 'eiv' ) THEN
582         IF( ktra == jp_tem )  htr_eiv(:,1) = ptr_sj( pva(:,:,:) )
583         IF( ktra == jp_sal )  str_eiv(:,1) = ptr_sj( pva(:,:,:) )
584      ENDIF
585      IF( cptr == 'vts' ) THEN
586         IF( ktra == jp_tem )  htr_vt(:,1) = ptr_sj( pva(:,:,:) )
587         IF( ktra == jp_sal )  str_vs(:,1) = ptr_sj( pva(:,:,:) )
588      ENDIF
589      !
590      IF( ln_subbas ) THEN
591         !
592         IF( cptr == 'adv' ) THEN
593             IF( ktra == jp_tem ) THEN
594                DO jn = 2, nptr
595                   htr_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
596                END DO
597             ENDIF
598             IF( ktra == jp_sal ) THEN
599                DO jn = 2, nptr
600                   str_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
601                END DO
602             ENDIF
603         ENDIF
604         IF( cptr == 'ldf' ) THEN
605             IF( ktra == jp_tem ) THEN
606                DO jn = 2, nptr
607                    htr_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
608                 END DO
609             ENDIF
610             IF( ktra == jp_sal ) THEN
611                DO jn = 2, nptr
612                   str_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
613                END DO
614             ENDIF
615         ENDIF
616         IF( cptr == 'eiv' ) THEN
617             IF( ktra == jp_tem ) THEN
618                DO jn = 2, nptr
619                    htr_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
620                 END DO
621             ENDIF
622             IF( ktra == jp_sal ) THEN
623                DO jn = 2, nptr
624                   str_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
625                END DO
626             ENDIF
627         ENDIF
628         IF( cptr == 'vts' ) THEN
629             IF( ktra == jp_tem ) THEN
630                DO jn = 2, nptr
631                    htr_vt(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
632                 END DO
633             ENDIF
634             IF( ktra == jp_sal ) THEN
635                DO jn = 2, nptr
636                   str_vs(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
637                END DO
638             ENDIF
639         ENDIF
640         !
641      ENDIF
642      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
643   END SUBROUTINE dia_ptr_ohst_components
644
645
646   FUNCTION dia_ptr_alloc()
647      !!----------------------------------------------------------------------
648      !!                    ***  ROUTINE dia_ptr_alloc  ***
649      !!----------------------------------------------------------------------
650      INTEGER               ::   dia_ptr_alloc   ! return value
651      INTEGER, DIMENSION(3) ::   ierr
652      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
653      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
654      REAL(KIND=jprb)               :: zhook_handle
655
656      CHARACTER(LEN=*), PARAMETER :: RoutineName='DIA_PTR_ALLOC'
657
658      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
659
660      !!----------------------------------------------------------------------
661      ierr(:) = 0
662      !
663      ALLOCATE( btmsk(jpi,jpj,nptr) ,              &
664         &      htr_adv(jpj,nptr) , str_adv(jpj,nptr) ,   &
665         &      htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) ,   &
666         &      htr_vt(jpj,nptr)  , str_vs(jpj,nptr)  ,   &
667         &      htr_ove(jpj,nptr) , str_ove(jpj,nptr) ,   &
668         &      htr_btr(jpj,nptr) , str_btr(jpj,nptr) ,   &
669         &      htr_ldf(jpj,nptr) , str_ldf(jpj,nptr) , STAT=ierr(1)  )
670         !
671      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2))
672      !
673      ALLOCATE( btm30(jpi,jpj), STAT=ierr(3)  )
674
675         !
676      dia_ptr_alloc = MAXVAL( ierr )
677      IF(lk_mpp)   CALL mpp_sum( dia_ptr_alloc )
678      !
679      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
680   END FUNCTION dia_ptr_alloc
681
682
683   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval )
684      !!----------------------------------------------------------------------
685      !!                    ***  ROUTINE ptr_sj_3d  ***
686      !!
687      !! ** Purpose :   i-k sum computation of a j-flux array
688      !!
689      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
690      !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
691      !!
692      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
693      !!----------------------------------------------------------------------
694      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pva   ! mask flux array at V-point
695      REAL(wp), INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
696      !
697      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
698      INTEGER                  ::   ijpj         ! ???
699      REAL(wp), POINTER, DIMENSION(:) :: p_fval  ! function value
700      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
701      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
702      REAL(KIND=jprb)               :: zhook_handle
703
704      CHARACTER(LEN=*), PARAMETER :: RoutineName='PTR_SJ_3D'
705
706      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
707
708      !!--------------------------------------------------------------------
709      !
710      p_fval => p_fval1d
711
712      ijpj = jpj
713      p_fval(:) = 0._wp
714      IF( PRESENT( pmsk ) ) THEN
715         DO jk = 1, jpkm1
716            DO jj = 2, jpjm1
717               DO ji = fs_2, fs_jpim1   ! Vector opt.
718                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) * pmsk(ji,jj)
719               END DO
720            END DO
721         END DO
722      ELSE
723         DO jk = 1, jpkm1
724            DO jj = 2, jpjm1
725               DO ji = fs_2, fs_jpim1   ! Vector opt.
726                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) 
727               END DO
728            END DO
729         END DO
730      ENDIF
731#if defined key_mpp_mpi
732      IF(lk_mpp)   CALL mpp_sum( p_fval, ijpj, ncomm_znl)
733#endif
734      !
735      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
736   END FUNCTION ptr_sj_3d
737
738
739   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval )
740      !!----------------------------------------------------------------------
741      !!                    ***  ROUTINE ptr_sj_2d  ***
742      !!
743      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array
744      !!
745      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
746      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
747      !!
748      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
749      !!----------------------------------------------------------------------
750      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pva   ! mask flux array at V-point
751      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
752      !
753      INTEGER                  ::   ji,jj       ! dummy loop arguments
754      INTEGER                  ::   ijpj        ! ???
755      REAL(wp), POINTER, DIMENSION(:) :: p_fval ! function value
756      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
757      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
758      REAL(KIND=jprb)               :: zhook_handle
759
760      CHARACTER(LEN=*), PARAMETER :: RoutineName='PTR_SJ_2D'
761
762      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
763
764      !!--------------------------------------------------------------------
765      !
766      p_fval => p_fval1d
767
768      ijpj = jpj
769      p_fval(:) = 0._wp
770      IF( PRESENT( pmsk ) ) THEN
771         DO jj = 2, jpjm1
772            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
773               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) * pmsk(ji,jj)
774            END DO
775         END DO
776      ELSE
777         DO jj = 2, jpjm1
778            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
779               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj)
780            END DO
781         END DO
782      ENDIF
783#if defined key_mpp_mpi
784      CALL mpp_sum( p_fval, ijpj, ncomm_znl )
785#endif
786      !
787      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
788   END FUNCTION ptr_sj_2d
789
790
791   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
792      !!----------------------------------------------------------------------
793      !!                    ***  ROUTINE ptr_sjk  ***
794      !!
795      !! ** Purpose :   i-sum computation of an array
796      !!
797      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
798      !!
799      !! ** Action  : - p_fval: i-mean poleward flux of pva
800      !!----------------------------------------------------------------------
801      !!
802      USE yomhook, ONLY: lhook, dr_hook
803      USE parkind1, ONLY: jprb, jpim
804
805      IMPLICIT none
806      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pta    ! mask flux array at V-point
807      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask
808      !!
809      INTEGER                           :: ji, jj, jk ! dummy loop arguments
810      REAL(wp), POINTER, DIMENSION(:,:) :: p_fval     ! return function value
811#if defined key_mpp_mpi
812      INTEGER, DIMENSION(1) ::   ish
813      INTEGER, DIMENSION(2) ::   ish2
814      INTEGER               ::   ijpjjpk
815      REAL(wp), DIMENSION(jpj*jpk) ::   zwork    ! mask flux array at V-point
816      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
817      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
818      REAL(KIND=jprb)               :: zhook_handle
819
820      CHARACTER(LEN=*), PARAMETER :: RoutineName='PTR_SJK'
821
822      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
823
824#endif
825      !!--------------------------------------------------------------------
826     !
827      p_fval => p_fval2d
828
829      p_fval(:,:) = 0._wp
830      !
831      IF( PRESENT( pmsk ) ) THEN
832         DO jk = 1, jpkm1
833            DO jj = 2, jpjm1
834!!gm here, use of tmask_i  ==> no need of loop over nldi, nlei....
835               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
836                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj)
837               END DO
838            END DO
839         END DO
840      ELSE
841         DO jk = 1, jpkm1
842            DO jj = 2, jpjm1
843               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
844                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * tmask_i(ji,jj)
845               END DO
846            END DO
847         END DO
848      END IF
849      !
850#if defined key_mpp_mpi
851      ijpjjpk = jpj*jpk
852      ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk
853      zwork(1:ijpjjpk) = RESHAPE( p_fval, ish )
854      CALL mpp_sum( zwork, ijpjjpk, ncomm_znl )
855      p_fval(:,:) = RESHAPE( zwork, ish2 )
856#endif
857      !
858
859      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
860   END FUNCTION ptr_sjk
861
862
863   !!======================================================================
864END MODULE diaptr
Note: See TracBrowser for help on using the repository browser.