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 branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

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

Last change on this file since 7179 was 7179, checked in by timgraham, 8 years ago

Manually merge in changes from v3.6_extra_CMIP6_diagnostics branch.
This change also includes a change of the domain_def.xml file so XIOS2 must be used from this revision onwards

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