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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90 @ 3255

Last change on this file since 3255 was 3255, checked in by charris, 12 years ago

Fixes for diaptr because the btm30 array was not being allocated, and the arrays for advective and diffusive heat and salt transport were not properly initialised. However there still appear to be unresolved problems with some of these diagnostics in my tests.

  • Property svn:keywords set to Id
File size: 43.1 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   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   dia_ptr      : Poleward Transport Diagnostics module
14   !!   dia_ptr_init : Initialization, namelist read
15   !!   dia_ptr_wri  : Output of poleward fluxes
16   !!   ptr_vjk      : "zonal" sum computation of a "meridional" flux array
17   !!   ptr_tjk      : "zonal" mean computation of a tracer field
18   !!   ptr_vj       : "zonal" and vertical sum computation of a "meridional" flux array
19   !!                   (Generic interface to ptr_vj_3d, ptr_vj_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       ! ocean active tracers: lateral physics
25   USE dianam           !
26   USE iom              ! IOM library
27   USE ioipsl           ! IO-IPSL library
28   USE in_out_manager   ! I/O manager
29   USE lib_mpp          ! MPP library
30   USE lbclnk           ! lateral boundary condition - processor exchanges
31   USE timing           ! preformance summary
32   USE wrk_nemo         ! working arrays
33
34   IMPLICIT NONE
35   PRIVATE
36
37   INTERFACE ptr_vj
38      MODULE PROCEDURE ptr_vj_3d, ptr_vj_2d
39   END INTERFACE
40
41   PUBLIC   dia_ptr_init   ! call in opa module
42   PUBLIC   dia_ptr        ! call in step module
43   PUBLIC   ptr_vj         ! call by tra_ldf & tra_adv routines
44   PUBLIC   ptr_vjk        ! call by tra_ldf & tra_adv routines
45
46   !                                           !!** namelist  namptr  **
47   LOGICAL , PUBLIC ::   ln_diaptr  = .FALSE.   !: Poleward transport flag (T) or not (F)
48   LOGICAL , PUBLIC ::   ln_subbas  = .FALSE.   !: Atlantic/Pacific/Indian basins calculation
49   LOGICAL , PUBLIC ::   ln_diaznl  = .FALSE.   !: Add zonal means and meridional stream functions
50   LOGICAL , PUBLIC ::   ln_ptrcomp = .FALSE.   !: Add decomposition : overturning (and gyre, soon ...)
51   INTEGER , PUBLIC ::   nn_fptr    = 15        !: frequency of ptr computation  [time step]
52   INTEGER , PUBLIC ::   nn_fwri    = 15        !: frequency of ptr outputs      [time step]
53
54   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf, htr_ove   !: Heat TRansports (adv, diff, overturn.)
55   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf, str_ove   !: Salt TRansports (adv, diff, overturn.)
56   
57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   btmsk                  ! T-point basin interior masks
58   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   btm30                  ! mask out Southern Ocean (=0 south of 30°S)
59   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htr  , str             ! adv heat and salt transports (approx)
60   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_jk, sn_jk , v_msf   ! i-mean T and S, j-Stream-Function
61   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sjk  , r1_sjk          ! i-mean i-k-surface and its inverse       
62   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htr_eiv, str_eiv       ! bolus adv heat ans salt transports ('key_diaeiv')
63   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_msf_eiv              ! bolus j-streamfuction              ('key_diaeiv')
64
65
66   INTEGER ::   niter       !
67   INTEGER ::   nidom_ptr   !
68   INTEGER ::   numptr      ! logical unit for Poleward TRansports
69   INTEGER ::   nptr        ! = 1 (ln_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (ln_subbas=T)
70
71   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
72   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rau0 x Cp)
73   REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg
74
75   REAL(wp), TARGET, DIMENSION(:),   ALLOCATABLE, SAVE :: p_fval1d
76   REAL(wp), TARGET, DIMENSION(:,:), ALLOCATABLE, SAVE :: p_fval2d
77
78   !! Integer, 1D workspace arrays. Not common enough to be implemented in
79   !! wrk_nemo module.
80   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex  , ndex_atl     , ndex_pac     , ndex_ind     , ndex_ipc
81   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::         ndex_atl_30  , ndex_pac_30  , ndex_ind_30  , ndex_ipc_30
82   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex_h, ndex_h_atl_30, ndex_h_pac_30, ndex_h_ind_30, ndex_h_ipc_30
83
84   !! * Substitutions
85#  include "domzgr_substitute.h90"
86#  include "vectopt_loop_substitute.h90"
87   !!----------------------------------------------------------------------
88   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
89   !! $Id$
90   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
91   !!----------------------------------------------------------------------
92CONTAINS
93
94   FUNCTION dia_ptr_alloc()
95      !!----------------------------------------------------------------------
96      !!                    ***  ROUTINE dia_ptr_alloc  ***
97      !!----------------------------------------------------------------------
98      INTEGER               ::   dia_ptr_alloc   ! return value
99      INTEGER, DIMENSION(6) ::   ierr
100      !!----------------------------------------------------------------------
101      ierr(:) = 0
102      !
103      ALLOCATE( btmsk(jpi,jpj,nptr) ,           &
104         &      htr_adv(jpj) , str_adv(jpj) ,   &
105         &      htr_ldf(jpj) , str_ldf(jpj) ,   &
106         &      htr_ove(jpj) , str_ove(jpj),    &
107         &      htr(jpj,nptr) , str(jpj,nptr) , &
108         &      tn_jk(jpj,jpk,nptr) , sn_jk (jpj,jpk,nptr) , v_msf(jpj,jpk,nptr) , &
109         &      sjk  (jpj,jpk,nptr) , r1_sjk(jpj,jpk,nptr) , STAT=ierr(1)  )
110         !
111#if defined key_diaeiv
112      ALLOCATE( htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) , &
113         &      v_msf_eiv(jpj,jpk,nptr) , STAT=ierr(2) )
114#endif
115      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(3))
116      !
117      ALLOCATE(ndex(jpj*jpk),        ndex_atl(jpj*jpk), ndex_pac(jpj*jpk), &
118         &     ndex_ind(jpj*jpk),    ndex_ipc(jpj*jpk),                    &
119         &     ndex_atl_30(jpj*jpk), ndex_pac_30(jpj*jpk), Stat=ierr(4))
120
121      ALLOCATE(ndex_ind_30(jpj*jpk), ndex_ipc_30(jpj*jpk),                   &
122         &     ndex_h(jpj),          ndex_h_atl_30(jpj), ndex_h_pac_30(jpj), &
123         &     ndex_h_ind_30(jpj),   ndex_h_ipc_30(jpj), Stat=ierr(5) )
124         !
125     ALLOCATE( btm30(jpi,jpj) , STAT=ierr(6)  )
126         !
127      dia_ptr_alloc = MAXVAL( ierr )
128      IF(lk_mpp)   CALL mpp_sum( dia_ptr_alloc )
129      !
130   END FUNCTION dia_ptr_alloc
131
132
133   FUNCTION ptr_vj_3d( pva )   RESULT ( p_fval )
134      !!----------------------------------------------------------------------
135      !!                    ***  ROUTINE ptr_vj_3d  ***
136      !!
137      !! ** Purpose :   i-k sum computation of a j-flux array
138      !!
139      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
140      !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
141      !!
142      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
143      !!----------------------------------------------------------------------
144      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pva   ! mask flux array at V-point
145      !!
146      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
147      INTEGER                  ::   ijpj         ! ???
148      REAL(wp), POINTER, DIMENSION(:) :: p_fval  ! function value
149      !!--------------------------------------------------------------------
150      !
151      p_fval => p_fval1d
152
153      ijpj = jpj
154      p_fval(:) = 0._wp
155      DO jk = 1, jpkm1
156         DO jj = 2, jpjm1
157            DO ji = fs_2, fs_jpim1   ! Vector opt.
158               p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) 
159            END DO
160         END DO
161      END DO
162#if defined key_mpp_mpi
163      IF(lk_mpp)   CALL mpp_sum( p_fval, ijpj, ncomm_znl)
164#endif
165      !
166   END FUNCTION ptr_vj_3d
167
168
169   FUNCTION ptr_vj_2d( pva )   RESULT ( p_fval )
170      !!----------------------------------------------------------------------
171      !!                    ***  ROUTINE ptr_vj_2d  ***
172      !!
173      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array
174      !!
175      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
176      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
177      !!
178      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
179      !!----------------------------------------------------------------------
180      IMPLICIT none
181      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point
182      !!
183      INTEGER                  ::   ji,jj       ! dummy loop arguments
184      INTEGER                  ::   ijpj        ! ???
185      REAL(wp), POINTER, DIMENSION(:) :: p_fval ! function value
186      !!--------------------------------------------------------------------
187      !
188      p_fval => p_fval1d
189
190      ijpj = jpj
191      p_fval(:) = 0._wp
192      DO jj = 2, jpjm1
193         DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
194            p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj)
195         END DO
196      END DO
197#if defined key_mpp_mpi
198      CALL mpp_sum( p_fval, ijpj, ncomm_znl )
199#endif
200      !
201   END FUNCTION ptr_vj_2d
202
203
204   FUNCTION ptr_vjk( pva, pmsk )   RESULT ( p_fval )
205      !!----------------------------------------------------------------------
206      !!                    ***  ROUTINE ptr_vjk  ***
207      !!
208      !! ** Purpose :   i-sum computation of a j-velocity array
209      !!
210      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
211      !!              pva is supposed to be a masked flux (i.e. * vmask)
212      !!
213      !! ** Action  : - p_fval: i-mean poleward flux of pva
214      !!----------------------------------------------------------------------
215      !!
216      IMPLICIT none
217      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pva    ! mask flux array at V-point
218      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask
219      !!
220      INTEGER                           :: ji, jj, jk ! dummy loop arguments
221      REAL(wp), POINTER, DIMENSION(:,:) :: p_fval     ! return function value
222#if defined key_mpp_mpi
223      INTEGER, DIMENSION(1) ::   ish
224      INTEGER, DIMENSION(2) ::   ish2
225      INTEGER               ::   ijpjjpk
226#endif
227#if defined key_mpp_mpi
228      REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point
229#endif
230      !!--------------------------------------------------------------------
231      !
232#if defined key_mpp_mpi
233      ijpjjpk = jpj*jpk
234      CALL wrk_alloc( jpj*jpk, zwork )
235#endif
236
237      p_fval => p_fval2d
238
239      p_fval(:,:) = 0._wp
240      !
241      IF( PRESENT( pmsk ) ) THEN
242         DO jk = 1, jpkm1
243            DO jj = 2, jpjm1
244!!gm here, use of tmask_i  ==> no need of loop over nldi, nlei....
245               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
246                  p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk) * pmsk(ji,jj)
247               END DO
248            END DO
249         END DO
250      ELSE
251         DO jk = 1, jpkm1
252            DO jj = 2, jpjm1
253               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
254                  p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji,jj)
255               END DO
256            END DO
257         END DO
258      END IF
259      !
260#if defined key_mpp_mpi
261      ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk
262      zwork(1:ijpjjpk) = RESHAPE( p_fval, ish )
263      CALL mpp_sum( zwork, ijpjjpk, ncomm_znl )
264      p_fval(:,:) = RESHAPE( zwork, ish2 )
265#endif
266      !
267#if defined key_mpp_mpi
268      CALL wrk_dealloc( jpj*jpk, zwork )
269#endif
270      !
271   END FUNCTION ptr_vjk
272
273
274   FUNCTION ptr_tjk( pta, pmsk )   RESULT ( p_fval )
275      !!----------------------------------------------------------------------
276      !!                    ***  ROUTINE ptr_tjk  ***
277      !!
278      !! ** Purpose :   i-sum computation of e1t*e3t * a tracer field
279      !!
280      !! ** Method  : - i-sum of mj(pta) using tmask
281      !!
282      !! ** Action  : - p_fval: i-sum of e1t*e3t*pta
283      !!----------------------------------------------------------------------
284      !!
285      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! tracer flux array at T-point
286      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)     ::   pmsk   ! Optional 2D basin mask
287      !!
288      INTEGER                           :: ji, jj, jk   ! dummy loop arguments
289      REAL(wp), POINTER, DIMENSION(:,:) :: p_fval       ! return function value
290#if defined key_mpp_mpi
291      INTEGER, DIMENSION(1) ::   ish
292      INTEGER, DIMENSION(2) ::   ish2
293      INTEGER               ::   ijpjjpk
294#endif
295#if defined key_mpp_mpi
296      REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point
297#endif
298      !!--------------------------------------------------------------------
299      !
300#if defined key_mpp_mpi
301      ijpjjpk = jpj*jpk
302      CALL wrk_alloc( jpj*jpk, zwork )
303#endif
304
305      p_fval => p_fval2d
306
307      p_fval(:,:) = 0._wp
308      DO jk = 1, jpkm1
309         DO jj = 2, jpjm1
310            DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
311               p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * e1t(ji,jj) * fse3t(ji,jj,jk) * pmsk(ji,jj)
312            END DO
313         END DO
314      END DO
315#if defined key_mpp_mpi
316      ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk
317      zwork(1:ijpjjpk)= RESHAPE( p_fval, ish )
318      CALL mpp_sum( zwork, ijpjjpk, ncomm_znl )
319      p_fval(:,:)= RESHAPE( zwork, ish2 )
320#endif
321      !
322#if defined key_mpp_mpi
323      CALL wrk_dealloc( jpj*jpk, zwork )
324#endif
325      !   
326   END FUNCTION ptr_tjk
327
328
329   SUBROUTINE dia_ptr( kt )
330      !!----------------------------------------------------------------------
331      !!                  ***  ROUTINE dia_ptr  ***
332      !!----------------------------------------------------------------------
333      USE oce,     vt  =>   ua   ! use ua as workspace
334      USE oce,     vs  =>   ua   ! use ua as workspace
335      IMPLICIT none
336      !!
337      INTEGER, INTENT(in) ::   kt   ! ocean time step index
338      !
339      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
340      REAL(wp) ::   zv               ! local scalar
341      !!----------------------------------------------------------------------
342      !
343      IF( nn_timing == 1 )   CALL timing_start('dia_ptr')
344      !
345      IF( kt == nit000 .OR. MOD( kt, nn_fptr ) == 0 )   THEN
346         !
347         IF( MOD( kt, nn_fptr ) == 0 ) THEN 
348            !
349            IF( ln_diaznl ) THEN               ! i-mean temperature and salinity
350               DO jn = 1, nptr
351                  tn_jk(:,:,jn) = ptr_tjk( tsn(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
352               END DO
353            ENDIF
354            !
355            !                          ! horizontal integral and vertical dz
356            !                                ! eulerian velocity
357            v_msf(:,:,1) = ptr_vjk( vn(:,:,:) ) 
358            DO jn = 2, nptr
359               v_msf(:,:,jn) = ptr_vjk( vn(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) 
360            END DO
361#if defined key_diaeiv
362            DO jn = 1, nptr                  ! bolus velocity
363               v_msf_eiv(:,:,jn) = ptr_vjk( v_eiv(:,:,:), btmsk(:,:,jn) )   ! here no btm30 for MSFeiv
364            END DO
365            !                                ! add bolus stream-function to the eulerian one
366            v_msf(:,:,:) = v_msf(:,:,:) + v_msf_eiv(:,:,:)
367#endif
368            !
369            !                          ! Transports
370            !                                ! local heat & salt transports at T-points  ( tsn*mj[vn+v_eiv] )
371            vt(:,:,jpk) = 0._wp   ;   vs(:,:,jpk) = 0._wp
372            DO jk= 1, jpkm1
373               DO jj = 2, jpj
374                  DO ji = 1, jpi
375#if defined key_diaeiv 
376                     zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) + v_eiv(ji,jj,jk) + v_eiv(ji,jj-1,jk) ) * 0.5_wp
377#else
378                     zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) ) * 0.5_wp
379#endif
380                     vt(:,jj,jk) = zv * tsn(:,jj,jk,jp_tem)
381                     vs(:,jj,jk) = zv * tsn(:,jj,jk,jp_sal)
382                  END DO
383               END DO
384            END DO
385!!gm useless as overlap areas are not used in ptr_vjk
386            CALL lbc_lnk( vs, 'V', -1. )   ;   CALL lbc_lnk( vt, 'V', -1. )
387!!gm
388            !                                ! heat & salt advective transports (approximation)
389            htr(:,1) = SUM( ptr_vjk( vt(:,:,:) ) , 2 ) * rc_pwatt   ! SUM over jk + conversion
390            str(:,1) = SUM( ptr_vjk( vs(:,:,:) ) , 2 ) * rc_ggram
391            DO jn = 2, nptr 
392               htr(:,jn) = SUM( ptr_vjk( vt(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) , 2 ) * rc_pwatt   ! mask Southern Ocean
393               str(:,jn) = SUM( ptr_vjk( vs(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) , 2 ) * rc_ggram   ! mask Southern Ocean
394            END DO
395
396            IF( ln_ptrcomp ) THEN            ! overturning transport
397               htr_ove(:) = SUM( v_msf(:,:,1) * tn_jk(:,:,1), 2 ) * rc_pwatt   ! SUM over jk + conversion
398               str_ove(:) = SUM( v_msf(:,:,1) * sn_jk(:,:,1), 2 ) * rc_ggram
399            END IF
400            !                                ! Advective and diffusive transport
401            htr_adv(:) = htr_adv(:) * rc_pwatt        ! these are computed in tra_adv... and tra_ldf... routines
402            htr_ldf(:) = htr_ldf(:) * rc_pwatt        ! here just the conversion in PW and Gg
403            str_adv(:) = str_adv(:) * rc_ggram
404            str_ldf(:) = str_ldf(:) * rc_ggram
405
406#if defined key_diaeiv
407            DO jn = 1, nptr                  ! Bolus component
408               htr_eiv(:,jn) = SUM( v_msf_eiv(:,:,jn) * tn_jk(:,:,jn), 2 ) * rc_pwatt   ! SUM over jk
409               str_eiv(:,jn) = SUM( v_msf_eiv(:,:,jn) * sn_jk(:,:,jn), 2 ) * rc_ggram   ! SUM over jk
410            END DO
411#endif
412            !                                ! "Meridional" Stream-Function
413            DO jn = 1, nptr
414               DO jk = 2, jpk 
415                  v_msf    (:,jk,jn) = v_msf    (:,jk-1,jn) + v_msf    (:,jk,jn)       ! Eulerian j-Stream-Function
416#if defined key_diaeiv
417                  v_msf_eiv(:,jk,jn) = v_msf_eiv(:,jk-1,jn) + v_msf_eiv(:,jk,jn)       ! Bolus    j-Stream-Function
418
419#endif
420               END DO
421            END DO
422            v_msf    (:,:,:) = v_msf    (:,:,:) * rc_sv       ! converte in Sverdrups
423#if defined key_diaeiv
424            v_msf_eiv(:,:,:) = v_msf_eiv(:,:,:) * rc_sv
425#endif
426         ENDIF
427         !
428         CALL dia_ptr_wri( kt )                        ! outputs
429         !
430      ENDIF
431      !
432#if defined key_mpp_mpi
433      IF( kt == nitend .AND. l_znl_root )   CALL histclo( numptr )      ! Close the file
434#else
435      IF( kt == nitend )                    CALL histclo( numptr )      ! Close the file
436#endif
437      !
438      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr')
439      !
440   END SUBROUTINE dia_ptr
441
442
443   SUBROUTINE dia_ptr_init
444      !!----------------------------------------------------------------------
445      !!                  ***  ROUTINE dia_ptr_init  ***
446      !!                   
447      !! ** Purpose :   Initialization, namelist read
448      !!----------------------------------------------------------------------
449      INTEGER ::   jn           ! dummy loop indices
450      INTEGER ::   inum, ierr   ! local integers
451#if defined key_mpp_mpi
452      INTEGER, DIMENSION(1) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
453#endif
454      !!
455      NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, ln_ptrcomp, nn_fptr, nn_fwri
456      !!----------------------------------------------------------------------
457      IF( nn_timing == 1 )   CALL timing_start('dia_ptr_init')
458
459      REWIND( numnam )                 ! Read Namelist namptr : poleward transport parameters
460      READ  ( numnam, 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,*) '      Overturning heat & salt transport                  ln_ptrcomp = ', ln_ptrcomp
469         WRITE(numout,*) '      T & S zonal mean and meridional stream function    ln_diaznl  = ', ln_diaznl 
470         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas
471         WRITE(numout,*) '      Frequency of computation                           nn_fptr    = ', nn_fptr
472         WRITE(numout,*) '      Frequency of outputs                               nn_fwri    = ', nn_fwri
473      ENDIF
474
475      IF( ln_subbas ) THEN   ;   nptr = 5       ! Global, Atlantic, Pacific, Indian, Indo-Pacific
476      ELSE                   ;   nptr = 1       ! Global only
477      ENDIF
478
479      !                                      ! allocate dia_ptr arrays
480      IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate arrays' )
481
482      rc_pwatt = rc_pwatt * rau0 * rcp          ! conversion from K.s-1 to PetaWatt
483
484      IF( .NOT. ln_diaptr ) THEN       ! diaptr not used
485        RETURN
486      ENDIF
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 )
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(:,:) = tmask(:,:,1)
499         END WHERE
500      ENDIF
501      btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
502     
503      DO jn = 1, nptr
504         btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
505      END DO
506     
507      IF( lk_vvl )   CALL ctl_stop( 'diaptr: error in vvl case as constant i-mean surface is used' )
508
509      !                                   ! i-sum of e1v*e3v surface and its inverse
510      DO jn = 1, nptr
511         sjk(:,:,jn) = ptr_tjk( tmask(:,:,:), btmsk(:,:,jn) )
512         r1_sjk(:,:,jn) = 0._wp
513         WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
514      END DO
515
516      ! Initialise arrays to zero because diatpr is called before they are first calculated
517      ! Note that this means diagnostics will not be exactly correct when model run is restarted.
518      htr_adv(:) = 0._wp ; str_adv(:) =  0._wp ;  htr_ldf(:) = 0._wp ; str_ldf(:) =  0._wp
519
520#if defined key_mpp_mpi 
521      iglo (1) = jpjglo                   ! MPP case using MPI  ('key_mpp_mpi')
522      iloc (1) = nlcj
523      iabsf(1) = njmppt(narea)
524      iabsl(:) = iabsf(:) + iloc(:) - 1
525      ihals(1) = nldj - 1
526      ihale(1) = nlcj - nlej
527      idid (1) = 2
528      CALL flio_dom_set( jpnj, nproc/jpni, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_ptr )
529#else
530      nidom_ptr = FLIO_DOM_NONE
531#endif
532      !
533      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr_init')
534      !
535   END SUBROUTINE dia_ptr_init
536
537
538   SUBROUTINE dia_ptr_wri( kt )
539      !!---------------------------------------------------------------------
540      !!                ***  ROUTINE dia_ptr_wri  ***
541      !!
542      !! ** Purpose :   output of poleward fluxes
543      !!
544      !! ** Method  :   NetCDF file
545      !!----------------------------------------------------------------------
546      !!
547      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
548      !!
549      INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw
550      INTEGER, SAVE ::   ndim  , ndim_atl     , ndim_pac     , ndim_ind     , ndim_ipc
551      INTEGER, SAVE ::           ndim_atl_30  , ndim_pac_30  , ndim_ind_30  , ndim_ipc_30
552      INTEGER, SAVE ::   ndim_h, ndim_h_atl_30, ndim_h_pac_30, ndim_h_ind_30, ndim_h_ipc_30
553      !!
554      CHARACTER (len=40) ::   clhstnam, clop, clop_once, cl_comment   ! temporary names
555      INTEGER            ::   iline, it, itmod, ji, jj, jk            !
556#if defined key_iomput
557      INTEGER            ::   inum                                    ! temporary logical unit
558#endif
559      REAL(wp)           ::   zsto, zout, zdt, zjulian                ! temporary scalars
560      !!
561      REAL(wp), POINTER, DIMENSION(:)   ::   zphi, zfoo    ! 1D workspace
562      REAL(wp), POINTER, DIMENSION(:,:) ::   z_1           ! 2D workspace
563      !!--------------------------------------------------------------------
564      !
565      CALL wrk_alloc( jpi      , zphi , zfoo )
566      CALL wrk_alloc( jpi , jpk, z_1 )
567
568      ! define time axis
569      it    = kt / nn_fptr
570      itmod = kt - nit000 + 1
571     
572      ! Initialization
573      ! --------------
574      IF( kt == nit000 ) THEN
575         niter = ( nit000 - 1 ) / nn_fptr
576         zdt = rdt
577         IF( nacc == 1 )   zdt = rdtmin
578         !
579         IF(lwp) THEN
580            WRITE(numout,*)
581            WRITE(numout,*) 'dia_ptr_wri : poleward transport and msf writing: initialization , niter = ', niter
582            WRITE(numout,*) '~~~~~~~~~~~~'
583         ENDIF
584
585         ! Reference latitude (used in plots)
586         ! ------------------
587         !                                           ! =======================
588         IF( cp_cfg == "orca" ) THEN                 !   ORCA configurations
589            !                                        ! =======================
590            IF( jp_cfg == 05  )   iline = 192   ! i-line that passes near the North Pole
591            IF( jp_cfg == 025 )   iline = 384   ! i-line that passes near the North Pole
592            IF( jp_cfg == 1   )   iline =  96   ! i-line that passes near the North Pole
593            IF( jp_cfg == 2   )   iline =  48   ! i-line that passes near the North Pole
594            IF( jp_cfg == 4   )   iline =  24   ! i-line that passes near the North Pole
595            zphi(1:jpj) = 0._wp
596            DO ji = mi0(iline), mi1(iline) 
597               zphi(1:jpj) = gphiv(ji,:)         ! if iline is in the local domain
598               ! Correct highest latitude for some configurations - will work if domain is parallelized in J ?
599               IF( jp_cfg == 05 ) THEN
600                  DO jj = mj0(jpjdta), mj1(jpjdta) 
601                     zphi( jj ) = zphi(mj0(jpjdta-1)) + ( zphi(mj0(jpjdta-1))-zphi(mj0(jpjdta-2)) ) * 0.5_wp
602                     zphi( jj ) = MIN( zphi(jj), 90._wp )
603                  END DO
604               END IF
605               IF( jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) THEN
606                  DO jj = mj0(jpjdta-1), mj1(jpjdta-1) 
607                     zphi( jj ) = 88.5_wp
608                  END DO
609                  DO jj = mj0(jpjdta  ), mj1(jpjdta  ) 
610                     zphi( jj ) = 89.5_wp
611                  END DO
612               END IF
613            END DO
614            ! provide the correct zphi to all local domains
615#if defined key_mpp_mpi
616            CALL mpp_sum( zphi, jpj, ncomm_znl )       
617#endif
618            !                                        ! =======================
619         ELSE                                        !   OTHER configurations
620            !                                        ! =======================
621            zphi(1:jpj) = gphiv(1,:)             ! assume lat/lon coordinate, select the first i-line
622            !
623         ENDIF
624         !
625         ! Work only on westmost processor (will not work if mppini2 is used)
626#if defined key_mpp_mpi
627         IF( l_znl_root ) THEN 
628#endif
629            !
630            ! OPEN netcdf file
631            ! ----------------
632            ! Define frequency of output and means
633            zsto = nn_fptr * zdt
634            IF( ln_mskland )   THEN    ! put 1.e+20 on land (very expensive!!)
635               clop      = "ave(only(x))"
636               clop_once = "once(only(x))"
637            ELSE                       ! no use of the mask value (require less cpu time)
638               clop      = "ave(x)"       
639               clop_once = "once"
640            ENDIF
641
642            zout = nn_fwri * zdt
643            zfoo(1:jpj) = 0._wp
644
645            CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )  ! Compute julian date from starting date of the run
646            zjulian = zjulian - adatrj                         ! set calendar origin to the beginning of the experiment
647
648#if defined key_iomput
649            ! Requested by IPSL people, use by their postpro...
650            IF(lwp) THEN
651               CALL dia_nam( clhstnam, nn_fwri,' ' )
652               CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
653               WRITE(inum,*) clhstnam
654               CLOSE(inum)
655            ENDIF
656#endif
657
658            CALL dia_nam( clhstnam, nn_fwri, 'diaptr' )
659            IF(lwp)WRITE( numout,*)" Name of diaptr NETCDF file : ", clhstnam
660
661            ! Horizontal grid : zphi()
662            CALL histbeg(clhstnam, 1, zfoo, jpj, zphi,   &
663               1, 1, 1, jpj, niter, zjulian, zdt*nn_fptr, nhoridz, numptr, domain_id=nidom_ptr)
664            ! Vertical grids : gdept_0, gdepw_0
665            CALL histvert( numptr, "deptht", "Vertical T levels",   &
666               &                   "m", jpk, gdept_0, ndepidzt, "down" )
667            CALL histvert( numptr, "depthw", "Vertical W levels",   &
668               &                   "m", jpk, gdepw_0, ndepidzw, "down" )
669            !
670            CALL wheneq ( jpj*jpk, MIN(sjk(:,:,1), 1._wp), 1, 1., ndex  , ndim  )      ! Lat-Depth
671            CALL wheneq ( jpj    , MIN(sjk(:,1,1), 1._wp), 1, 1., ndex_h, ndim_h )     ! Lat
672
673            IF( ln_subbas ) THEN
674               z_1(:,1) = 1._wp
675               WHERE ( gphit(jpi/2,:) < -30._wp )   z_1(:,1) = 0._wp
676               DO jk = 2, jpk
677                  z_1(:,jk) = z_1(:,1)
678               END DO
679               !                       ! Atlantic (jn=2)
680               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,2)         , 1._wp), 1, 1., ndex_atl     , ndim_atl      ) ! Lat-Depth
681               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,2)*z_1(:,:), 1._wp), 1, 1., ndex_atl_30  , ndim_atl_30   ) ! Lat-Depth
682               CALL wheneq ( jpj    , MIN(sjk(:,1,2)*z_1(:,1), 1._wp), 1, 1., ndex_h_atl_30, ndim_h_atl_30 ) ! Lat
683               !                       ! Pacific (jn=3)
684               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,3)         , 1._wp), 1, 1., ndex_pac     , ndim_pac      ) ! Lat-Depth
685               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,3)*z_1(:,:), 1._wp), 1, 1., ndex_pac_30  , ndim_pac_30   ) ! Lat-Depth
686               CALL wheneq ( jpj    , MIN(sjk(:,1,3)*z_1(:,1), 1._wp), 1, 1., ndex_h_pac_30, ndim_h_pac_30 ) ! Lat
687               !                       ! Indian (jn=4)
688               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,4)         , 1._wp), 1, 1., ndex_ind     , ndim_ind      ) ! Lat-Depth
689               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,4)*z_1(:,:), 1._wp), 1, 1., ndex_ind_30  , ndim_ind_30   ) ! Lat-Depth
690               CALL wheneq ( jpj    , MIN(sjk(:,1,4)*z_1(:,1), 1._wp), 1, 1., ndex_h_ind_30, ndim_h_ind_30 ) ! Lat
691               !                       ! Indo-Pacific (jn=5)
692               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,5)         , 1._wp), 1, 1., ndex_ipc     , ndim_ipc      ) ! Lat-Depth
693               CALL wheneq ( jpj*jpk, MIN(sjk(:,:,5)*z_1(:,:), 1._wp), 1, 1., ndex_ipc_30  , ndim_ipc_30   ) ! Lat-Depth
694               CALL wheneq ( jpj    , MIN(sjk(:,1,5)*z_1(:,1), 1._wp), 1, 1., ndex_h_ipc_30, ndim_h_ipc_30 ) ! Lat
695            ENDIF
696            !
697#if defined key_diaeiv
698            cl_comment = ' (Bolus part included)'
699#else
700            cl_comment = '                      '
701#endif
702            IF( ln_diaznl ) THEN             !  Zonal mean T and S
703               CALL histdef( numptr, "zotemglo", "Zonal Mean Temperature","C" ,   &
704                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
705               CALL histdef( numptr, "zosalglo", "Zonal Mean Salinity","PSU"  ,   &
706                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
707
708               CALL histdef( numptr, "zosrfglo", "Zonal Mean Surface","m^2"   ,   &
709                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
710               !
711               IF (ln_subbas) THEN
712                  CALL histdef( numptr, "zotematl", "Zonal Mean Temperature: Atlantic","C" ,   &
713                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
714                  CALL histdef( numptr, "zosalatl", "Zonal Mean Salinity: Atlantic","PSU"  ,   &
715                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
716                  CALL histdef( numptr, "zosrfatl", "Zonal Mean Surface: Atlantic","m^2"   ,   &
717                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
718
719                  CALL histdef( numptr, "zotempac", "Zonal Mean Temperature: Pacific","C"  ,   &
720                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
721                  CALL histdef( numptr, "zosalpac", "Zonal Mean Salinity: Pacific","PSU"   ,   &
722                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
723                  CALL histdef( numptr, "zosrfpac", "Zonal Mean Surface: Pacific","m^2"    ,   &
724                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
725
726                  CALL histdef( numptr, "zotemind", "Zonal Mean Temperature: Indian","C"   ,   &
727                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
728                  CALL histdef( numptr, "zosalind", "Zonal Mean Salinity: Indian","PSU"    ,   &
729                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
730                  CALL histdef( numptr, "zosrfind", "Zonal Mean Surface: Indian","m^2"     ,   &
731                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
732
733                  CALL histdef( numptr, "zotemipc", "Zonal Mean Temperature: Pacific+Indian","C" ,   &
734                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
735                  CALL histdef( numptr, "zosalipc", "Zonal Mean Salinity: Pacific+Indian","PSU"  ,   &
736                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
737                  CALL histdef( numptr, "zosrfipc", "Zonal Mean Surface: Pacific+Indian","m^2"   ,   &
738                     1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
739               ENDIF
740            ENDIF
741            !
742            !  Meridional Stream-Function (Eulerian and Bolus)
743            CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: Global"//TRIM(cl_comment),"Sv" ,   &
744               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
745            IF( ln_subbas .AND. ln_diaznl ) THEN
746               CALL histdef( numptr, "zomsfatl", "Meridional Stream-Function: Atlantic"//TRIM(cl_comment),"Sv" ,   &
747                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
748               CALL histdef( numptr, "zomsfpac", "Meridional Stream-Function: Pacific"//TRIM(cl_comment),"Sv"  ,   &
749                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
750               CALL histdef( numptr, "zomsfind", "Meridional Stream-Function: Indian"//TRIM(cl_comment),"Sv"   ,   &
751                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
752               CALL histdef( numptr, "zomsfipc", "Meridional Stream-Function: Indo-Pacific"//TRIM(cl_comment),"Sv" ,&
753                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
754            ENDIF
755            !
756            !  Heat transport
757            CALL histdef( numptr, "sophtadv", "Advective Heat Transport"      ,   &
758               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
759            CALL histdef( numptr, "sophtldf", "Diffusive Heat Transport"      ,   &
760               "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
761            IF ( ln_ptrcomp ) THEN
762               CALL histdef( numptr, "sophtove", "Overturning Heat Transport"    ,   &
763                  "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
764            END IF
765            IF( ln_subbas ) THEN
766               CALL histdef( numptr, "sohtatl", "Heat Transport Atlantic"//TRIM(cl_comment),  &
767                  "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
768               CALL histdef( numptr, "sohtpac", "Heat Transport Pacific"//TRIM(cl_comment) ,  &
769                  "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
770               CALL histdef( numptr, "sohtind", "Heat Transport Indian"//TRIM(cl_comment)  ,  &
771                  "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
772               CALL histdef( numptr, "sohtipc", "Heat Transport Pacific+Indian"//TRIM(cl_comment), &
773                  "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
774            ENDIF
775            !
776            !  Salt transport
777            CALL histdef( numptr, "sopstadv", "Advective Salt Transport"      ,   &
778               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
779            CALL histdef( numptr, "sopstldf", "Diffusive Salt Transport"      ,   &
780               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
781            IF ( ln_ptrcomp ) THEN
782               CALL histdef( numptr, "sopstove", "Overturning Salt Transport"    ,   &
783                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
784            END IF
785#if defined key_diaeiv
786            ! Eddy induced velocity
787            CALL histdef( numptr, "zomsfeiv", "Bolus Meridional Stream-Function: global",   &
788               "Sv"      , 1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
789            CALL histdef( numptr, "sophteiv", "Bolus Advective Heat Transport",   &
790               "PW"      , 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
791            CALL histdef( numptr, "sopsteiv", "Bolus Advective Salt Transport",   &
792               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
793#endif
794            IF( ln_subbas ) THEN
795               CALL histdef( numptr, "sostatl", "Salt Transport Atlantic"//TRIM(cl_comment)      ,  &
796                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
797               CALL histdef( numptr, "sostpac", "Salt Transport Pacific"//TRIM(cl_comment)      ,   &
798                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
799               CALL histdef( numptr, "sostind", "Salt Transport Indian"//TRIM(cl_comment)      ,    &
800                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
801               CALL histdef( numptr, "sostipc", "Salt Transport Pacific+Indian"//TRIM(cl_comment),  &
802                  "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
803            ENDIF
804            !
805            CALL histend( numptr )
806            !
807         END IF
808#if defined key_mpp_mpi
809      END IF
810#endif
811
812#if defined key_mpp_mpi
813      IF( MOD( itmod, nn_fptr ) == 0 .AND. l_znl_root ) THEN
814#else
815      IF( MOD( itmod, nn_fptr ) == 0  ) THEN
816#endif
817         niter = niter + 1
818
819         IF( ln_diaznl ) THEN
820            CALL histwrite( numptr, "zosrfglo", niter, sjk  (:,:,1) , ndim, ndex )
821            CALL histwrite( numptr, "zotemglo", niter, tn_jk(:,:,1)  , ndim, ndex )
822            CALL histwrite( numptr, "zosalglo", niter, sn_jk(:,:,1)  , ndim, ndex )
823
824            IF (ln_subbas) THEN
825               CALL histwrite( numptr, "zosrfatl", niter, sjk(:,:,2), ndim_atl, ndex_atl )
826               CALL histwrite( numptr, "zosrfpac", niter, sjk(:,:,3), ndim_pac, ndex_pac )
827               CALL histwrite( numptr, "zosrfind", niter, sjk(:,:,4), ndim_ind, ndex_ind )
828               CALL histwrite( numptr, "zosrfipc", niter, sjk(:,:,5), ndim_ipc, ndex_ipc )
829
830               CALL histwrite( numptr, "zotematl", niter, tn_jk(:,:,2)  , ndim_atl, ndex_atl )
831               CALL histwrite( numptr, "zosalatl", niter, sn_jk(:,:,2)  , ndim_atl, ndex_atl )
832               CALL histwrite( numptr, "zotempac", niter, tn_jk(:,:,3)  , ndim_pac, ndex_pac )
833               CALL histwrite( numptr, "zosalpac", niter, sn_jk(:,:,3)  , ndim_pac, ndex_pac )
834               CALL histwrite( numptr, "zotemind", niter, tn_jk(:,:,4)  , ndim_ind, ndex_ind )
835               CALL histwrite( numptr, "zosalind", niter, sn_jk(:,:,4)  , ndim_ind, ndex_ind )
836               CALL histwrite( numptr, "zotemipc", niter, tn_jk(:,:,5)  , ndim_ipc, ndex_ipc )
837               CALL histwrite( numptr, "zosalipc", niter, sn_jk(:,:,5)  , ndim_ipc, ndex_ipc )
838            END IF
839         ENDIF
840
841         ! overturning outputs:
842         CALL histwrite( numptr, "zomsfglo", niter, v_msf(:,:,1), ndim, ndex )
843         IF( ln_subbas .AND. ln_diaznl ) THEN
844            CALL histwrite( numptr, "zomsfatl", niter, v_msf(:,:,2) , ndim_atl_30, ndex_atl_30 )
845            CALL histwrite( numptr, "zomsfpac", niter, v_msf(:,:,3) , ndim_pac_30, ndex_pac_30 )
846            CALL histwrite( numptr, "zomsfind", niter, v_msf(:,:,4) , ndim_ind_30, ndex_ind_30 )
847            CALL histwrite( numptr, "zomsfipc", niter, v_msf(:,:,5) , ndim_ipc_30, ndex_ipc_30 )
848         ENDIF
849#if defined key_diaeiv
850         CALL histwrite( numptr, "zomsfeiv", niter, v_msf_eiv(:,:,1), ndim  , ndex   )
851#endif
852
853         ! heat transport outputs:
854         IF( ln_subbas ) THEN
855            CALL histwrite( numptr, "sohtatl", niter, htr(:,2)  , ndim_h_atl_30, ndex_h_atl_30 )
856            CALL histwrite( numptr, "sohtpac", niter, htr(:,3)  , ndim_h_pac_30, ndex_h_pac_30 )
857            CALL histwrite( numptr, "sohtind", niter, htr(:,4)  , ndim_h_ind_30, ndex_h_ind_30 )
858            CALL histwrite( numptr, "sohtipc", niter, htr(:,5)  , ndim_h_ipc_30, ndex_h_ipc_30 )
859            CALL histwrite( numptr, "sostatl", niter, str(:,2)  , ndim_h_atl_30, ndex_h_atl_30 )
860            CALL histwrite( numptr, "sostpac", niter, str(:,3)  , ndim_h_pac_30, ndex_h_pac_30 )
861            CALL histwrite( numptr, "sostind", niter, str(:,4)  , ndim_h_ind_30, ndex_h_ind_30 )
862            CALL histwrite( numptr, "sostipc", niter, str(:,5)  , ndim_h_ipc_30, ndex_h_ipc_30 )
863         ENDIF
864
865         CALL histwrite( numptr, "sophtadv", niter, htr_adv     , ndim_h, ndex_h )
866         CALL histwrite( numptr, "sophtldf", niter, htr_ldf     , ndim_h, ndex_h )
867         CALL histwrite( numptr, "sopstadv", niter, str_adv     , ndim_h, ndex_h )
868         CALL histwrite( numptr, "sopstldf", niter, str_ldf     , ndim_h, ndex_h )
869         IF( ln_ptrcomp ) THEN
870            CALL histwrite( numptr, "sopstove", niter, str_ove(:) , ndim_h, ndex_h )
871            CALL histwrite( numptr, "sophtove", niter, htr_ove(:) , ndim_h, ndex_h )
872         ENDIF
873#if defined key_diaeiv
874         CALL histwrite( numptr, "sophteiv", niter, htr_eiv(:,1)  , ndim_h, ndex_h )
875         CALL histwrite( numptr, "sopsteiv", niter, str_eiv(:,1)  , ndim_h, ndex_h )
876#endif
877         !
878      ENDIF
879      !
880      CALL wrk_dealloc( jpi      , zphi , zfoo )
881      CALL wrk_dealloc( jpi , jpk, z_1 )
882      !
883  END SUBROUTINE dia_ptr_wri
884
885   !!======================================================================
886END MODULE diaptr
Note: See TracBrowser for help on using the repository browser.