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

source: branches/UKMO/dev_r5518_GO6_package_diaptr_gm_fix/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90 @ 8606

Last change on this file since 8606 was 8606, checked in by timgraham, 7 years ago

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