source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90 @ 8638

Last change on this file since 8638 was 8638, checked in by frrh, 3 years ago

Applied fix from Tim Graham's branh: branches/UKMO/dev_r5518_GO6_package_msfeiv_fix
at revision 8635.

Met Office GMED ticket 356 refers.

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