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 NEMO/trunk/src/OCE/DIA – NEMO

source: NEMO/trunk/src/OCE/DIA/diaptr.F90 @ 14994

Last change on this file since 14994 was 14834, checked in by hadcv, 3 years ago

#2600: Merge in dev_r14273_HPC-02_Daley_Tiling

  • Property svn:keywords set to Id
File size: 34.1 KB
RevLine 
[134]1MODULE diaptr
2   !!======================================================================
3   !!                       ***  MODULE  diaptr  ***
[1340]4   !! Ocean physics:  Computes meridonal transports and zonal means
[134]5   !!=====================================================================
[1559]6   !! History :  1.0  ! 2003-09  (C. Talandier, G. Madec)  Original code
7   !!            2.0  ! 2006-01  (A. Biastoch)  Allow sub-basins computation
[2528]8   !!            3.2  ! 2010-03  (O. Marti, S. Flavoni) Add fields
9   !!            3.3  ! 2010-10  (G. Madec)  dynamical allocation
[5147]10   !!            3.6  ! 2014-12  (C. Ethe) use of IOM
[7646]11   !!            3.6  ! 2016-06  (T. Graham) Addition of diagnostics for CMIP6
[12276]12   !!            4.0  ! 2010-08  ( C. Ethe, J. Deshayes ) Improvment
[134]13   !!----------------------------------------------------------------------
[508]14
15   !!----------------------------------------------------------------------
[134]16   !!   dia_ptr      : Poleward Transport Diagnostics module
17   !!   dia_ptr_init : Initialization, namelist read
[5147]18   !!   ptr_sjk      : "zonal" mean computation of a field - tracer or flux array
19   !!   ptr_sj       : "zonal" and vertical sum computation of a "meridional" flux array
20   !!                   (Generic interface to ptr_sj_3d, ptr_sj_2d)
[134]21   !!----------------------------------------------------------------------
[2528]22   USE oce              ! ocean dynamics and active tracers
23   USE dom_oce          ! ocean space and time domain
[14090]24   USE domtile
[2528]25   USE phycst           ! physical constants
[5147]26   !
[2528]27   USE iom              ! IOM library
28   USE in_out_manager   ! I/O manager
29   USE lib_mpp          ! MPP library
[3294]30   USE timing           ! preformance summary
[134]31
32   IMPLICIT NONE
33   PRIVATE
34
[13982]35   INTERFACE ptr_sum
36      MODULE PROCEDURE ptr_sum_3d, ptr_sum_2d
37   END INTERFACE
38
[5147]39   INTERFACE ptr_sj
40      MODULE PROCEDURE ptr_sj_3d, ptr_sj_2d
[134]41   END INTERFACE
42
[508]43   PUBLIC   dia_ptr        ! call in step module
[7646]44   PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines
[134]45
[13982]46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hstr_adv, hstr_ldf, hstr_eiv   !: Heat/Salt TRansports(adv, diff, Bolus.)
47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional)
48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   pvtr_int, pzon_int             !: Other zonal integrals
[1345]49
[13982]50   LOGICAL, PUBLIC    ::   l_diaptr       !: tracers  trend flag
51   INTEGER, PARAMETER ::   jp_msk = 3
52   INTEGER, PARAMETER ::   jp_vtr = 4
[2715]53
[2528]54   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
[12489]55   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rho0 x Cp)
56   REAL(wp) ::   rc_ggram = 1.e-9_wp   ! conversion from g    to Gg  (further x rho0)
[134]57
[12276]58   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
59   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk34 ! mask out Southern Ocean (=0 south of 34°S)
[2715]60
[13982]61   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag
[2715]62
[134]63   !! * Substitutions
[12377]64#  include "do_loop_substitute.h90"
[13237]65#  include "domzgr_substitute.h90"
[134]66   !!----------------------------------------------------------------------
[9598]67   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[14072]68   !! $Id$
[10068]69   !! Software governed by the CeCILL license (see ./LICENSE)
[134]70   !!----------------------------------------------------------------------
71CONTAINS
72
[14834]73   ! NOTE: [tiling] tiling sometimes changes the diagnostics very slightly, usually where there are few zonal points e.g. the northern Indian Ocean basin. The difference is usually very small, for one point in one diagnostic. Presumably this is because of the additional zonal integration step over tiles.
[12377]74   SUBROUTINE dia_ptr( kt, Kmm, pvtr )
[5147]75      !!----------------------------------------------------------------------
76      !!                  ***  ROUTINE dia_ptr  ***
77      !!----------------------------------------------------------------------
[14072]78      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index
[12377]79      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
[13982]80      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
81      !!----------------------------------------------------------------------
[5147]82      !
[13982]83      IF( ln_timing )   CALL timing_start('dia_ptr')
84
85      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init    ! -> will define l_diaptr and nbasin
86      !
87      IF( l_diaptr ) THEN
88         ! Calculate zonal integrals
89         IF( PRESENT( pvtr ) ) THEN
90            CALL dia_ptr_zint( Kmm, pvtr )
91         ELSE
92            CALL dia_ptr_zint( Kmm )
93         ENDIF
94
95         ! Calculate diagnostics only when zonal integrals have finished
[14834]96         IF( .NOT. l_istiled .OR. ntile == nijtile ) CALL dia_ptr_iom(kt, Kmm, pvtr)
[13982]97      ENDIF
98
99      IF( ln_timing )   CALL timing_stop('dia_ptr')
100      !
101   END SUBROUTINE dia_ptr
102
103
104   SUBROUTINE dia_ptr_iom( kt, Kmm, pvtr )
105      !!----------------------------------------------------------------------
106      !!                  ***  ROUTINE dia_ptr_iom  ***
107      !!----------------------------------------------------------------------
108      !! ** Purpose : Calculate diagnostics and send to XIOS
109      !!----------------------------------------------------------------------
110      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index
111      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
112      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
113      !
[5147]114      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
115      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
[12276]116      REAL(wp), DIMENSION(jpj)      ::  zvsum, ztsum, zssum   ! 1D workspace
[7646]117      !
118      !overturning calculation
[13557]119      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse
120      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   zt_jk, zs_jk        ! i-mean T and S, j-Stream-Function
[7646]121
[13557]122      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   z4d1, z4d2
123      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   z3dtr
[5147]124      !!----------------------------------------------------------------------
125      !
[13982]126      ALLOCATE( z3dtr(jpi,jpj,nbasin) )
[12377]127
[5147]128      IF( PRESENT( pvtr ) ) THEN
[12276]129         IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF
[13557]130            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin) )
[13982]131            !
[13557]132            DO jn = 1, nbasin                                    ! by sub-basins
[13982]133               z4d1(1,:,:,jn) =  pvtr_int(:,:,jp_vtr,jn)                  ! zonal cumulative effective transport excluding closed seas
134               DO jk = jpkm1, 1, -1
[12276]135                  z4d1(1,:,jk,jn) = z4d1(1,:,jk+1,jn) - z4d1(1,:,jk,jn)    ! effective j-Stream-Function (MSF)
[5147]136               END DO
[13982]137               DO ji = 2, jpi
[12276]138                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
[5147]139               ENDDO
140            END DO
[12276]141            CALL iom_put( 'zomsf', z4d1 * rc_sv )
[13982]142            !
[13557]143            DEALLOCATE( z4d1 )
[5147]144         ENDIF
[12276]145         IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
[13982]146            ALLOCATE( sjk(jpj,jpk,nbasin), r1_sjk(jpj,jpk,nbasin), v_msf(jpj,jpk,nbasin),   &
147               &      zt_jk(jpj,jpk,nbasin), zs_jk(jpj,jpk,nbasin) )
148            !
[13557]149            DO jn = 1, nbasin
[13982]150               sjk(:,:,jn) = pvtr_int(:,:,jp_msk,jn)
[12276]151               r1_sjk(:,:,jn) = 0._wp
152               WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
153               ! i-mean T and S, j-Stream-Function, basin
[13982]154               zt_jk(:,:,jn) = pvtr_int(:,:,jp_tem,jn) * r1_sjk(:,:,jn)
155               zs_jk(:,:,jn) = pvtr_int(:,:,jp_sal,jn) * r1_sjk(:,:,jn)
156               v_msf(:,:,jn) = pvtr_int(:,:,jp_vtr,jn)
[12276]157               hstr_ove(:,jp_tem,jn) = SUM( v_msf(:,:,jn)*zt_jk(:,:,jn), 2 )
158               hstr_ove(:,jp_sal,jn) = SUM( v_msf(:,:,jn)*zs_jk(:,:,jn), 2 )
159               !
160            ENDDO
[13557]161            DO jn = 1, nbasin
[12276]162               z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
[13982]163               DO ji = 2, jpi
[12276]164                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
165               ENDDO
166            ENDDO
167            CALL iom_put( 'sophtove', z3dtr )
[13557]168            DO jn = 1, nbasin
[12276]169               z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
[13982]170               DO ji = 2, jpi
[12276]171                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
172               ENDDO
173            ENDDO
174            CALL iom_put( 'sopstove', z3dtr )
[13982]175            !
176            DEALLOCATE( sjk, r1_sjk, v_msf, zt_jk, zs_jk )
[12276]177         ENDIF
[11993]178
[12276]179         IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
[14072]180            ! Calculate barotropic heat and salt transport here
[13982]181            ALLOCATE( sjk(jpj,1,nbasin), r1_sjk(jpj,1,nbasin) )
182            !
[13557]183            DO jn = 1, nbasin
[13982]184               sjk(:,1,jn) = SUM( pvtr_int(:,:,jp_msk,jn), 2 )
[12276]185               r1_sjk(:,1,jn) = 0._wp
186               WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
187               !
[13982]188               zvsum(:) =    SUM( pvtr_int(:,:,jp_vtr,jn), 2 )
189               ztsum(:) =    SUM( pvtr_int(:,:,jp_tem,jn), 2 )
190               zssum(:) =    SUM( pvtr_int(:,:,jp_sal,jn), 2 )
[12276]191               hstr_btr(:,jp_tem,jn) = zvsum(:) * ztsum(:) * r1_sjk(:,1,jn)
192               hstr_btr(:,jp_sal,jn) = zvsum(:) * zssum(:) * r1_sjk(:,1,jn)
193               !
[7646]194            ENDDO
[13557]195            DO jn = 1, nbasin
[12276]196               z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
[13982]197               DO ji = 2, jpi
[12276]198                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
199               ENDDO
[7646]200            ENDDO
[12276]201            CALL iom_put( 'sophtbtr', z3dtr )
[13557]202            DO jn = 1, nbasin
[12276]203               z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
[13982]204               DO ji = 2, jpi
[12276]205                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
[7646]206               ENDDO
[12276]207            ENDDO
208            CALL iom_put( 'sopstbtr', z3dtr )
[13982]209            !
210            DEALLOCATE( sjk, r1_sjk )
211         ENDIF
[5147]212         !
[13982]213         hstr_ove(:,:,:) = 0._wp       ! Zero before next timestep
214         hstr_btr(:,:,:) = 0._wp
215         pvtr_int(:,:,:,:) = 0._wp
[5147]216      ELSE
[13982]217         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface
[13557]218            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin), z4d2(jpi,jpj,jpk,nbasin) )
[12276]219            !
[13557]220            DO jn = 1, nbasin
[13982]221               z4d1(1,:,:,jn) = pzon_int(:,:,jp_msk,jn)
222               DO ji = 2, jpi
223                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
[13557]224               ENDDO
[12276]225            ENDDO
226            CALL iom_put( 'zosrf', z4d1 )
227            !
[13557]228            DO jn = 1, nbasin
[13982]229               z4d2(1,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
230               DO ji = 2, jpi
[12276]231                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
[5147]232               ENDDO
[12276]233            ENDDO
234            CALL iom_put( 'zotem', z4d2 )
235            !
[13557]236            DO jn = 1, nbasin
[13982]237               z4d2(1,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
238               DO ji = 2, jpi
[12276]239                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
[5147]240               ENDDO
[12276]241            ENDDO
242            CALL iom_put( 'zosal', z4d2 )
[13982]243            !
[13557]244            DEALLOCATE( z4d1, z4d2 )
[5147]245         ENDIF
246         !
247         !                                ! Advective and diffusive heat and salt transport
[14072]248         IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN
249            !
[13557]250            DO jn = 1, nbasin
[12276]251               z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
[13982]252               DO ji = 2, jpi
[12276]253                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
254               ENDDO
[11993]255            ENDDO
[12276]256            CALL iom_put( 'sophtadv', z3dtr )
[13557]257            DO jn = 1, nbasin
[12276]258               z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
[13982]259               DO ji = 2, jpi
[12276]260                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
261               ENDDO
[11993]262            ENDDO
[12276]263            CALL iom_put( 'sopstadv', z3dtr )
264         ENDIF
265         !
[14072]266         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN
267            !
[13557]268            DO jn = 1, nbasin
[12276]269               z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
[13982]270               DO ji = 2, jpi
[12276]271                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
[11989]272               ENDDO
[12276]273            ENDDO
274            CALL iom_put( 'sophtldf', z3dtr )
[13557]275            DO jn = 1, nbasin
[12276]276               z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
[13982]277               DO ji = 2, jpi
[12276]278                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
[11989]279               ENDDO
[12276]280            ENDDO
281            CALL iom_put( 'sopstldf', z3dtr )
[11989]282         ENDIF
283         !
[14072]284         IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN
285            !
[13557]286            DO jn = 1, nbasin
[12276]287               z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
[13982]288               DO ji = 2, jpi
[12276]289                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
[7646]290               ENDDO
[12276]291            ENDDO
292            CALL iom_put( 'sophteiv', z3dtr )
[13557]293            DO jn = 1, nbasin
[12276]294               z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
[13982]295               DO ji = 2, jpi
[12276]296                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
[7646]297               ENDDO
[12276]298            ENDDO
299            CALL iom_put( 'sopsteiv', z3dtr )
[11993]300         ENDIF
[12276]301         !
302         IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
[13557]303             DO jn = 1, nbasin
[12276]304                z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
[13982]305                DO ji = 2, jpi
[12276]306                   z3dtr(ji,:,jn) = z3dtr(1,:,jn)
307                ENDDO
308             ENDDO
309             CALL iom_put( 'sophtvtr', z3dtr )
[13557]310             DO jn = 1, nbasin
[12276]311               z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
[13982]312               DO ji = 2, jpi
[12276]313                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
314               ENDDO
315            ENDDO
316            CALL iom_put( 'sopstvtr', z3dtr )
[7646]317         ENDIF
[5147]318         !
[12276]319         IF( iom_use( 'uocetr_vsum_cumul' ) ) THEN
320            CALL iom_get_var(  'uocetr_vsum_op', z2d ) ! get uocetr_vsum_op from xml
[14072]321            z2d(:,:) = ptr_ci_2d( z2d(:,:) )
[12276]322            CALL iom_put( 'uocetr_vsum_cumul', z2d )
323         ENDIF
324         !
[13982]325         hstr_adv(:,:,:) = 0._wp       ! Zero before next timestep
326         hstr_ldf(:,:,:) = 0._wp
327         hstr_eiv(:,:,:) = 0._wp
328         hstr_vtr(:,:,:) = 0._wp
329         pzon_int(:,:,:,:) = 0._wp
[5147]330      ENDIF
331      !
[13557]332      DEALLOCATE( z3dtr )
333      !
[13982]334   END SUBROUTINE dia_ptr_iom
[5147]335
336
[13982]337   SUBROUTINE dia_ptr_zint( Kmm, pvtr )
338      !!----------------------------------------------------------------------
339      !!                    ***  ROUTINE dia_ptr_zint ***
340      !!----------------------------------------------------------------------
341      !! ** Purpose : i and i-k sum operations on arrays
342      !!
343      !! ** Method  : - Call ptr_sjk (i sum) or ptr_sj (i-k sum) to perform the sum operation
344      !!              - Call ptr_sum to add this result to the sum over tiles
345      !!
346      !! ** Action  : pvtr_int - terms for volume streamfunction, heat/salt transport barotropic/overturning terms
347      !!              pzon_int - terms for i mean temperature/salinity
348      !!----------------------------------------------------------------------
349      INTEGER                     , INTENT(in)           :: Kmm          ! time level index
350      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in), OPTIONAL :: pvtr         ! j-effective transport
351      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: zmask        ! 3D workspace
352      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          :: zts          ! 4D workspace
353      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: sjk, v_msf   ! Zonal sum: i-k surface area, j-effective transport
354      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: zt_jk, zs_jk ! Zonal sum: i-k surface area * (T, S)
355      REAL(wp)                                           :: zsfc, zvfc   ! i-k surface area
356      INTEGER  ::   ji, jj, jk, jn                                       ! dummy loop indices
357      !!----------------------------------------------------------------------
358
359      IF( PRESENT( pvtr ) ) THEN
360         ! i sum of effective j transport excluding closed seas
361         IF( iom_use( 'zomsf' ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
362            ALLOCATE( v_msf(A1Dj(nn_hls),jpk,nbasin) )
363
364            DO jn = 1, nbasin
365               v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )
366            ENDDO
367
368            CALL ptr_sum( pvtr_int(:,:,jp_vtr,:), v_msf(:,:,:) )
369
370            DEALLOCATE( v_msf )
371         ENDIF
372
373         ! i sum of j surface area, j surface area - temperature/salinity product on V grid
374         IF(  iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.   &
375            & iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
376            ALLOCATE( zmask(A2D(nn_hls),jpk), zts(A2D(nn_hls),jpk,jpts), &
377               &      sjk(A1Dj(nn_hls),jpk,nbasin), &
378               &      zt_jk(A1Dj(nn_hls),jpk,nbasin), zs_jk(A1Dj(nn_hls),jpk,nbasin) )
379
380            zmask(:,:,:) = 0._wp
381            zts(:,:,:,:) = 0._wp
382
[14215]383            DO_3D( 1, 1, 1, 0, 1, jpkm1 )
[13982]384               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
385               zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
386               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc !Tracers averaged onto V grid
387               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
388            END_3D
389
390            DO jn = 1, nbasin
391               sjk(:,:,jn)   = ptr_sjk( zmask(:,:,:)     , btmsk(:,:,jn) )
392               zt_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
393               zs_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
394            ENDDO
395
396            CALL ptr_sum( pvtr_int(:,:,jp_msk,:), sjk(:,:,:)   )
397            CALL ptr_sum( pvtr_int(:,:,jp_tem,:), zt_jk(:,:,:) )
398            CALL ptr_sum( pvtr_int(:,:,jp_sal,:), zs_jk(:,:,:) )
399
400            DEALLOCATE( zmask, zts, sjk, zt_jk, zs_jk )
401         ENDIF
402      ELSE
403         ! i sum of j surface area - temperature/salinity product on T grid
404         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN
405            ALLOCATE( zmask(A2D(nn_hls),jpk), zts(A2D(nn_hls),jpk,jpts), &
406               &      sjk(A1Dj(nn_hls),jpk,nbasin), &
407               &      zt_jk(A1Dj(nn_hls),jpk,nbasin), zs_jk(A1Dj(nn_hls),jpk,nbasin) )
408
409            zmask(:,:,:) = 0._wp
410            zts(:,:,:,:) = 0._wp
411
412            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
413               zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm)
414               zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
415               zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc
416               zts(ji,jj,jk,jp_sal) = ts(ji,jj,jk,jp_sal,Kmm) * zsfc
417            END_3D
418
419            DO jn = 1, nbasin
420               sjk(:,:,jn)   = ptr_sjk( zmask(:,:,:)     , btmsk(:,:,jn) )
421               zt_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
422               zs_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
423            ENDDO
424
425            CALL ptr_sum( pzon_int(:,:,jp_msk,:), sjk(:,:,:)   )
426            CALL ptr_sum( pzon_int(:,:,jp_tem,:), zt_jk(:,:,:) )
427            CALL ptr_sum( pzon_int(:,:,jp_sal,:), zs_jk(:,:,:) )
428
429            DEALLOCATE( zmask, zts, sjk, zt_jk, zs_jk )
430         ENDIF
431
432         ! i-k sum of j surface area - temperature/salinity product on V grid
433         IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
434            ALLOCATE( zts(A2D(nn_hls),jpk,jpts) )
435
436            zts(:,:,:,:) = 0._wp
437
[14215]438            DO_3D( 1, 1, 1, 0, 1, jpkm1 )
[13982]439               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
440               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid
441               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
442            END_3D
443
444            CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) )
445            CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) )
446
447            DEALLOCATE( zts )
448         ENDIF
449      ENDIF
450   END SUBROUTINE dia_ptr_zint
451
452
[5147]453   SUBROUTINE dia_ptr_init
454      !!----------------------------------------------------------------------
455      !!                  ***  ROUTINE dia_ptr_init  ***
[14072]456      !!
[13557]457      !! ** Purpose :   Initialization
[5147]458      !!----------------------------------------------------------------------
[12377]459      INTEGER ::  inum, jn           ! local integers
[5147]460      !!
[12276]461      REAL(wp), DIMENSION(jpi,jpj) :: zmsk
[5147]462      !!----------------------------------------------------------------------
[13982]463
[13557]464      ! l_diaptr is defined with iom_use
465      !   --> dia_ptr_init must be done after the call to iom_init
466      !   --> cannot be .TRUE. without cpp key: key_iom -->  nbasin define by iom_init is initialized
467      l_diaptr = iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  &
468         &       iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  &
469         &       iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  &
[13982]470         &       iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  &
[13557]471         &       iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  &
[13982]472         &       iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' )
[14072]473
[5147]474      IF(lwp) THEN                     ! Control print
475         WRITE(numout,*)
476         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
477         WRITE(numout,*) '~~~~~~~~~~~~'
[12377]478         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr
[5147]479      ENDIF
480
[14072]481      IF( l_diaptr ) THEN
[5147]482         !
483         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
[13557]484         !
[12489]485         rc_pwatt = rc_pwatt * rho0_rcp          ! conversion from K.s-1 to PetaWatt
486         rc_ggram = rc_ggram * rho0              ! conversion from m3/s to Gg/s
[5147]487
488         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
489
[14072]490         btmsk(:,:,1) = tmask_i(:,:)
[13557]491         IF( nbasin == 5 ) THEN   ! nbasin has been initialized in iom_init to define the axis "basin"
492            CALL iom_open( 'subbasins', inum )
493            CALL iom_get( inum, jpdom_global, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
494            CALL iom_get( inum, jpdom_global, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
495            CALL iom_get( inum, jpdom_global, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
496            CALL iom_close( inum )
497            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )            ! Indo-Pacific basin
498         ENDIF
499         DO jn = 2, nbasin
500            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)                 ! interior domain only
[5147]501         END DO
[12276]502         ! JD : modification so that overturning streamfunction is available in Atlantic at 34S to compare with observations
503         WHERE( gphit(:,:)*tmask_i(:,:) < -34._wp)
504           zmsk(:,:) = 0._wp      ! mask out Southern Ocean
[14072]505         ELSE WHERE
[12276]506           zmsk(:,:) = ssmask(:,:)
507         END WHERE
[14072]508         btmsk34(:,:,1) = btmsk(:,:,1)
[13557]509         DO jn = 2, nbasin
510            btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:)                  ! interior domain only
[12276]511         ENDDO
[5147]512
513         ! Initialise arrays to zero because diatpr is called before they are first calculated
514         ! Note that this means diagnostics will not be exactly correct when model run is restarted.
[14072]515         hstr_adv(:,:,:) = 0._wp
516         hstr_ldf(:,:,:) = 0._wp
517         hstr_eiv(:,:,:) = 0._wp
518         hstr_ove(:,:,:) = 0._wp
[12276]519         hstr_btr(:,:,:) = 0._wp           !
520         hstr_vtr(:,:,:) = 0._wp           !
[13982]521         pvtr_int(:,:,:,:) = 0._wp
522         pzon_int(:,:,:,:) = 0._wp
[7753]523         !
[12377]524         ll_init = .FALSE.
525         !
[14072]526      ENDIF
527      !
[5147]528   END SUBROUTINE dia_ptr_init
529
[9124]530
[14072]531   SUBROUTINE dia_ptr_hst( ktra, cptr, pvflx )
[7646]532      !!----------------------------------------------------------------------
533      !!                    ***  ROUTINE dia_ptr_hst ***
534      !!----------------------------------------------------------------------
535      !! Wrapper for heat and salt transport calculations to calculate them for each basin
536      !! Called from all advection and/or diffusion routines
537      !!----------------------------------------------------------------------
538      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
539      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
[13982]540      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: pvflx ! 3D input array of advection/diffusion
541      REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin)                 :: zsj   !
[7646]542      INTEGER                                        :: jn    !
[5147]543
[13982]544      DO jn = 1, nbasin
545         zsj(:,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
546      ENDDO
[12276]547      !
[7646]548      IF( cptr == 'adv' ) THEN
[13982]549         IF( ktra == jp_tem )  CALL ptr_sum( hstr_adv(:,jp_tem,:), zsj(:,:) )
550         IF( ktra == jp_sal )  CALL ptr_sum( hstr_adv(:,jp_sal,:), zsj(:,:) )
551      ELSE IF( cptr == 'ldf' ) THEN
552         IF( ktra == jp_tem )  CALL ptr_sum( hstr_ldf(:,jp_tem,:), zsj(:,:) )
553         IF( ktra == jp_sal )  CALL ptr_sum( hstr_ldf(:,jp_sal,:), zsj(:,:) )
554      ELSE IF( cptr == 'eiv' ) THEN
555         IF( ktra == jp_tem )  CALL ptr_sum( hstr_eiv(:,jp_tem,:), zsj(:,:) )
556         IF( ktra == jp_sal )  CALL ptr_sum( hstr_eiv(:,jp_sal,:), zsj(:,:) )
557      ELSE IF( cptr == 'vtr' ) THEN
558         IF( ktra == jp_tem )  CALL ptr_sum( hstr_vtr(:,jp_tem,:), zsj(:,:) )
559         IF( ktra == jp_sal )  CALL ptr_sum( hstr_vtr(:,jp_sal,:), zsj(:,:) )
[7646]560      ENDIF
[12276]561      !
[13982]562   END SUBROUTINE dia_ptr_hst
563
564
565   SUBROUTINE ptr_sum_2d( phstr, pva )
566      !!----------------------------------------------------------------------
567      !!                    ***  ROUTINE ptr_sum_2d ***
568      !!----------------------------------------------------------------------
569      !! ** Purpose : Add two 2D arrays with (j,nbasin) dimensions
570      !!
571      !! ** Method  : - phstr = phstr + pva
572      !!              - Call mpp_sum if the final tile
573      !!
574      !! ** Action  : phstr
575      !!----------------------------------------------------------------------
576      REAL(wp), DIMENSION(jpj,nbasin) , INTENT(inout)         ::  phstr  !
577      REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin), INTENT(in)            ::  pva    !
578      INTEGER                                               ::  jj
[14229]579#if ! defined key_mpi_off
[13982]580      INTEGER, DIMENSION(1)           ::  ish1d
581      INTEGER, DIMENSION(2)           ::  ish2d
582      REAL(wp), DIMENSION(jpj*nbasin) ::  zwork
583#endif
584
585      DO jj = ntsj, ntej
586         phstr(jj,:) = phstr(jj,:)  + pva(jj,:)
587      END DO
588
[14229]589#if ! defined key_mpi_off
[14834]590      IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
[13982]591         ish1d(1) = jpj*nbasin
592         ish2d(1) = jpj ; ish2d(2) = nbasin
593         zwork(:) = RESHAPE( phstr(:,:), ish1d )
594         CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
595         phstr(:,:) = RESHAPE( zwork, ish2d )
[7646]596      ENDIF
[13982]597#endif
598   END SUBROUTINE ptr_sum_2d
599
600
601   SUBROUTINE ptr_sum_3d( phstr, pva )
602      !!----------------------------------------------------------------------
603      !!                    ***  ROUTINE ptr_sum_3d ***
604      !!----------------------------------------------------------------------
605      !! ** Purpose : Add two 3D arrays with (j,k,nbasin) dimensions
606      !!
607      !! ** Method  : - phstr = phstr + pva
608      !!              - Call mpp_sum if the final tile
609      !!
610      !! ** Action  : phstr
611      !!----------------------------------------------------------------------
612      REAL(wp), DIMENSION(jpj,jpk,nbasin) , INTENT(inout)     ::  phstr  !
613      REAL(wp), DIMENSION(A1Dj(nn_hls),jpk,nbasin), INTENT(in)        ::  pva    !
614      INTEGER                                               ::  jj, jk
[14229]615#if ! defined key_mpi_off
[13982]616      INTEGER, DIMENSION(1)              ::  ish1d
617      INTEGER, DIMENSION(3)              ::  ish3d
618      REAL(wp), DIMENSION(jpj*jpk*nbasin)  ::  zwork
619#endif
620
621      DO jk = 1, jpk
622         DO jj = ntsj, ntej
623            phstr(jj,jk,:) = phstr(jj,jk,:)  + pva(jj,jk,:)
624         END DO
625      END DO
626
[14229]627#if ! defined key_mpi_off
[14834]628      IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
[13982]629         ish1d(1) = jpj*jpk*nbasin
630         ish3d(1) = jpj ; ish3d(2) = jpk ; ish3d(3) = nbasin
631         zwork(:) = RESHAPE( phstr(:,:,:), ish1d )
632         CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
633         phstr(:,:,:) = RESHAPE( zwork, ish3d )
[7646]634      ENDIF
[13982]635#endif
636   END SUBROUTINE ptr_sum_3d
[7646]637
638
[2715]639   FUNCTION dia_ptr_alloc()
640      !!----------------------------------------------------------------------
641      !!                    ***  ROUTINE dia_ptr_alloc  ***
642      !!----------------------------------------------------------------------
643      INTEGER               ::   dia_ptr_alloc   ! return value
[13982]644      INTEGER, DIMENSION(2) ::   ierr
[2715]645      !!----------------------------------------------------------------------
646      ierr(:) = 0
647      !
[13557]648      ! nbasin has been initialized in iom_init to define the axis "basin"
649      !
[12276]650      IF( .NOT. ALLOCATED( btmsk ) ) THEN
[13557]651         ALLOCATE( btmsk(jpi,jpj,nbasin)    , btmsk34(jpi,jpj,nbasin),   &
652            &      hstr_adv(jpj,jpts,nbasin), hstr_eiv(jpj,jpts,nbasin), &
653            &      hstr_ove(jpj,jpts,nbasin), hstr_btr(jpj,jpts,nbasin), &
654            &      hstr_ldf(jpj,jpts,nbasin), hstr_vtr(jpj,jpts,nbasin), STAT=ierr(1)  )
[12276]655            !
[13982]656         ALLOCATE( pvtr_int(jpj,jpk,jpts+2,nbasin), &
657            &      pzon_int(jpj,jpk,jpts+1,nbasin), STAT=ierr(2) )
[2715]658         !
[12276]659         dia_ptr_alloc = MAXVAL( ierr )
660         CALL mpp_sum( 'diaptr', dia_ptr_alloc )
661      ENDIF
[2715]662      !
663   END FUNCTION dia_ptr_alloc
664
665
[12377]666   FUNCTION ptr_sj_3d( pvflx, pmsk )   RESULT ( p_fval )
[134]667      !!----------------------------------------------------------------------
[5147]668      !!                    ***  ROUTINE ptr_sj_3d  ***
[134]669      !!
[2528]670      !! ** Purpose :   i-k sum computation of a j-flux array
[134]671      !!
[12377]672      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
673      !!              pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
[134]674      !!
[12377]675      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
[508]676      !!----------------------------------------------------------------------
[13982]677      REAL(wp), INTENT(in), DIMENSION(A2D(nn_hls),jpk)  ::   pvflx  ! mask flux array at V-point
678      REAL(wp), INTENT(in), DIMENSION(jpi,jpj)  ::   pmsk   ! Optional 2D basin mask
[5147]679      !
[508]680      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
[13982]681      REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval  ! function value
[134]682      !!--------------------------------------------------------------------
[508]683      !
[2528]684      p_fval(:) = 0._wp
[13295]685      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]686         p_fval(jj) = p_fval(jj) + pvflx(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
687      END_3D
[5147]688   END FUNCTION ptr_sj_3d
[134]689
690
[12377]691   FUNCTION ptr_sj_2d( pvflx, pmsk )   RESULT ( p_fval )
[134]692      !!----------------------------------------------------------------------
[5147]693      !!                    ***  ROUTINE ptr_sj_2d  ***
[134]694      !!
[12377]695      !! ** Purpose :   "zonal" and vertical sum computation of a j-flux array
[134]696      !!
[12377]697      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
698      !!      pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
[134]699      !!
[12377]700      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
[508]701      !!----------------------------------------------------------------------
[13982]702      REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls))     ::   pvflx  ! mask flux array at V-point
[12276]703      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
[5147]704      !
[2715]705      INTEGER                  ::   ji,jj       ! dummy loop arguments
[13982]706      REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval ! function value
[134]707      !!--------------------------------------------------------------------
[13982]708      !
[2528]709      p_fval(:) = 0._wp
[13295]710      DO_2D( 0, 0, 0, 0 )
[12377]711         p_fval(jj) = p_fval(jj) + pvflx(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj)
712      END_2D
[5147]713   END FUNCTION ptr_sj_2d
[134]714
[12276]715   FUNCTION ptr_ci_2d( pva )   RESULT ( p_fval )
716      !!----------------------------------------------------------------------
717      !!                    ***  ROUTINE ptr_ci_2d  ***
718      !!
719      !! ** Purpose :   "meridional" cumulated sum computation of a j-flux array
720      !!
721      !! ** Method  : - j cumulated sum of pva using the interior 2D vmask (umask_i).
722      !!
723      !! ** Action  : - p_fval: j-cumulated sum of pva
724      !!----------------------------------------------------------------------
725      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)  ::   pva   ! mask flux array at V-point
726      !
727      INTEGER                  ::   ji,jj,jc       ! dummy loop arguments
[14072]728      INTEGER                  ::   ijpj        ! ???
[12276]729      REAL(wp), DIMENSION(jpi,jpj) :: p_fval ! function value
730      !!--------------------------------------------------------------------
[14072]731      !
[12276]732      ijpj = jpj  ! ???
733      p_fval(:,:) = 0._wp
734      DO jc = 1, jpnj ! looping over all processors in j axis
[13295]735         DO_2D( 0, 0, 0, 0 )
[12377]736            p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj)
737         END_2D
[12276]738      END DO
[14072]739      !
[12276]740   END FUNCTION ptr_ci_2d
[134]741
[12276]742
743
[5147]744   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
[134]745      !!----------------------------------------------------------------------
[5147]746      !!                    ***  ROUTINE ptr_sjk  ***
[134]747      !!
[5147]748      !! ** Purpose :   i-sum computation of an array
[134]749      !!
[12377]750      !! ** Method  : - i-sum of field using the interior 2D vmask (pmsk).
[134]751      !!
[12377]752      !! ** Action  : - p_fval: i-sum of masked field
[508]753      !!----------------------------------------------------------------------
[2715]754      !!
755      IMPLICIT none
[13982]756      REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls),jpk) ::   pta    ! mask flux array at V-point
757      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
[134]758      !!
[2715]759      INTEGER                           :: ji, jj, jk ! dummy loop arguments
[13982]760      REAL(wp), DIMENSION(A1Dj(nn_hls),jpk) :: p_fval     ! return function value
[134]761      !!--------------------------------------------------------------------
[1559]762      !
[2528]763      p_fval(:,:) = 0._wp
[508]764      !
[13295]765      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]766         p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
767      END_3D
[5147]768   END FUNCTION ptr_sjk
[134]769
[1559]770
[134]771   !!======================================================================
772END MODULE diaptr
Note: See TracBrowser for help on using the repository browser.