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/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

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