Ticket #361: diaptr.F90

File diaptr.F90, 40.2 KB (added by marti, 12 years ago)

diaptr.F90

Line 
1MODULE diaptr
2   !!======================================================================
3   !!                       ***  MODULE  diaptr  ***
4   !! Ocean physics:  Computes meridonal transports and zonal means
5   !!=====================================================================
6   !! History :  9.0  !  03-09  (C. Talandier, G. Madec)  Original code
7   !!            9.0  !  06-01  (A. Biastoch)  Allow sub-basins computation
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   dia_ptr      : Poleward Transport Diagnostics module
12   !!   dia_ptr_init : Initialization, namelist read
13   !!   dia_ptr_wri  : Output of poleward fluxes
14   !!   ptr_vjk      : "zonal" sum computation of a "meridional" flux array
15   !!   ptr_tjk      : "zonal" mean computation of a tracer field
16   !!   ptr_vj       : "zonal" and vertical sum computation of a "meridional"
17   !!                : flux array; Generic interface: ptr_vj_3d, ptr_vj_2d
18   !!----------------------------------------------------------------------
19   USE oce           ! ocean dynamics and active tracers
20   USE dom_oce       ! ocean space and time domain
21   USE ldftra_oce    ! ???
22   USE lib_mpp
23   USE in_out_manager
24   USE dianam
25   USE phycst
26   USE iom
27   USE ioipsl         
28   USE daymod
29
30   IMPLICIT NONE
31   PRIVATE
32
33   INTERFACE ptr_vj
34      MODULE PROCEDURE ptr_vj_3d, ptr_vj_2d
35   END INTERFACE
36
37   PUBLIC   dia_ptr_init   ! call in opa module
38   PUBLIC   dia_ptr        ! call in step module
39   PUBLIC   ptr_vj         ! call by tra_ldf & tra_adv routines
40   PUBLIC   ptr_vjk        ! call by tra_ldf & tra_adv routines
41
42!!! ** init namelist (namptr)
43   LOGICAL , PUBLIC                 ::   ln_diaptr = .FALSE.   !: Poleward transport flag (T) or not (F)
44   LOGICAL , PUBLIC                 ::   ln_subbas = .FALSE.   !: Atlantic/Pacific/Indian basins calculation
45   LOGICAL , PUBLIC                 ::   ln_diaznl             !: Add zonal means and meridional stream functions
46   INTEGER , PUBLIC                 ::   nf_ptr = 15           !: frequency of ptr computation
47   INTEGER , PUBLIC                 ::   nf_ptr_wri = 15       !: frequency of ptr outputs
48
49   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_adv, pst_adv      !: heat and salt poleward transport: advection
50   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_ove, pst_ove      !: heat and salt poleward transport: overturning
51   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_ldf, pst_ldf      !: heat and salt poleward transport: lateral diffusion
52#if defined key_diaeiv
53   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_eiv, pst_eiv      !: heat and salt poleward transport: bolus advection
54#endif
55   REAL(wp), PUBLIC, DIMENSION(jpj) ::   ht_glo,ht_atl,ht_ind,ht_pac,ht_ipc !: heat
56   REAL(wp), PUBLIC, DIMENSION(jpj) ::   st_glo,st_atl,st_ind,st_pac,st_ipc !: salt
57
58   REAL(wp), DIMENSION(jpj,jpk) ::   tn_jk_glo, sn_jk_glo,  &  !: "zonal" mean temperature and salinity
59      &                              tn_jk_atl, sn_jk_atl,  &
60      &                              tn_jk_pac, sn_jk_pac,  &
61      &                              tn_jk_ind, sn_jk_ind,  &
62      &                              tn_jk_ipc, sn_jk_ipc,  &
63      &                              v_msf_glo       ,  &  !: "meridional" Stream-Function
64      &                              v_msf_atl       ,  & 
65      &                              v_msf_pac       ,  & 
66      &                              v_msf_ind       ,  & 
67      &                              v_msf_ipc       ,  & 
68      &                              surf_jk_glo     ,  &  !: Ocean "zonal" section surface
69      &                              surf_jk_atl     ,  & 
70      &                              surf_jk_pac     ,  & 
71      &                              surf_jk_ind     ,  &         
72      &                              surf_jk_ipc     ,  & 
73      &                              surf_jk_r_glo   ,  &  !: inverse of the ocean "zonal" section surface
74      &                              surf_jk_r_atl   ,  & 
75      &                              surf_jk_r_pac   ,  & 
76      &                              surf_jk_r_ind   ,  &         
77      &                              surf_jk_r_ipc       
78#if defined key_diaeiv
79   REAL(wp), DIMENSION(jpj,jpk) ::   v_msf_eiv              !: bolus "meridional" Stream-Function
80#endif
81   REAL(wp), DIMENSION(jpi,jpj) ::   abasin, pbasin, ibasin, dbasin, sbasin !: Sub basin masks
82
83   !! * Substitutions
84#  include "domzgr_substitute.h90"
85#  include "vectopt_loop_substitute.h90"
86   !!----------------------------------------------------------------------
87   !!   OPA 9.0 , LOCEAN-IPSL (2005)
88   !! $Id: diaptr.F90 1317 2009-02-17 16:31:42Z smasson $
89   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
90   !!----------------------------------------------------------------------
91
92CONTAINS
93
94   FUNCTION ptr_vj_3d( pva )   RESULT ( p_fval )
95      !!----------------------------------------------------------------------
96      !!                    ***  ROUTINE ptr_vj_3d  ***
97      !!
98      !! ** Purpose :   "zonal" and vertical sum computation of a "meridional"
99      !!      flux array
100      !!
101      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
102      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
103      !!
104      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
105      !!----------------------------------------------------------------------
106      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pva   ! mask flux array at V-point
107      !!
108      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
109      INTEGER                  ::   ijpj         ! ???
110      REAL(wp), DIMENSION(jpj) ::   p_fval       ! function value
111      !!--------------------------------------------------------------------
112      !
113      ijpj = jpj
114      p_fval(:) = 0.e0
115      DO jk = 1, jpkm1
116         DO jj = 2, jpjm1
117            DO ji = fs_2, fs_jpim1   ! Vector opt.
118               p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) 
119            END DO
120         END DO
121      END DO
122      !
123      IF( lk_mpp )   CALL mpp_sum( p_fval, ijpj )     !!bug  I presume
124      !
125   END FUNCTION ptr_vj_3d
126
127
128   FUNCTION ptr_vj_2d( pva )   RESULT ( p_fval )
129      !!----------------------------------------------------------------------
130      !!                    ***  ROUTINE ptr_vj_2d  ***
131      !!
132      !! ** Purpose :   "zonal" and vertical sum computation of a "meridional"
133      !!      flux array
134      !!
135      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
136      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
137      !!
138      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
139      !!----------------------------------------------------------------------
140      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point
141      !!
142      INTEGER                  ::   ji,jj    ! dummy loop arguments
143      INTEGER                  ::   ijpj     ! ???
144      REAL(wp), DIMENSION(jpj) ::   p_fval   ! function value
145      !!--------------------------------------------------------------------
146      !
147      ijpj = jpj
148      p_fval(:) = 0.e0
149      DO jj = 2, jpjm1
150         DO ji = fs_2, fs_jpim1   ! Vector opt.
151            p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj+1) * tmask_i(ji,jj)
152         END DO
153      END DO
154      !
155      IF( lk_mpp )   CALL mpp_sum( p_fval, ijpj )     !!bug  I presume
156      !
157   END FUNCTION ptr_vj_2d
158
159
160   FUNCTION ptr_vjk( pva, bmask )   RESULT ( p_fval )
161      !!----------------------------------------------------------------------
162      !!                    ***  ROUTINE ptr_vjk  ***
163      !!
164      !! ** Purpose :   "zonal" sum computation of a "meridional" flux array
165      !!
166      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
167      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
168      !!
169      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
170      !!----------------------------------------------------------------------
171      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pva   ! mask flux array at V-point
172      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL :: bmask ! Optional 2D basin mask
173      !!
174      INTEGER                      ::   ji, jj, jk   ! dummy loop arguments
175      INTEGER , DIMENSION (1)      ::   ish
176      INTEGER , DIMENSION (2)      ::   ish2
177      REAL(wp), DIMENSION(jpj*jpk) ::   zwork        ! temporary vector for mpp_sum
178      REAL(wp), DIMENSION(jpj,jpk) ::   p_fval       ! return function value
179      !!--------------------------------------------------------------------
180      !
181      p_fval(:,:) = 0.e0
182      !
183      IF (PRESENT (bmask)) THEN
184         DO jk = 1, jpkm1
185            DO jj = 2, jpjm1
186               DO ji = fs_2, fs_jpim1
187                  p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk)   &
188                     &                                          * tmask_i(ji,jj+1) * tmask_i(ji,jj) &
189                     &                                          * bmask(ji,jj)
190               END DO
191            END DO
192         END DO
193      ELSE
194         DO jk = 1, jpkm1
195            DO jj = 2, jpjm1
196               DO ji = fs_2, fs_jpim1
197                  p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk)   &
198                     &                                          * tmask_i(ji,jj+1) * tmask_i(ji,jj)
199               END DO
200            END DO
201         END DO
202      END IF
203      !
204      IF(lk_mpp) THEN
205         ish(1) = jpj*jpk  ;  ish2(1) = jpj  ;  ish2(2) = jpk
206         zwork(:)= RESHAPE( p_fval, ish )
207         CALL mpp_sum( zwork, jpj*jpk )
208         p_fval(:,:)= RESHAPE( zwork, ish2 )
209      END IF
210      !
211   END FUNCTION ptr_vjk
212
213   FUNCTION ptr_tjk( pta, bmask )   RESULT ( p_fval )
214      !!----------------------------------------------------------------------
215      !!                    ***  ROUTINE ptr_tjk  ***
216      !!
217      !! ** Purpose :   "zonal" mean computation of a tracer field
218      !!
219      !! ** Method  : - i-sum of mj(pta) using tmask
220      !!      multiplied by the inverse of the surface of the "zonal" ocean
221      !!      section
222      !!
223      !! ** Action  : - p_fval: i-k-mean poleward flux of pta
224      !!----------------------------------------------------------------------
225      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! tracer flux array at T-point
226      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL :: bmask ! Optional 2D basin mask
227      !!
228      INTEGER                     ::   ji, jj, jk   ! dummy loop arguments
229      INTEGER, DIMENSION (1)      ::   ish
230      INTEGER, DIMENSION (2)      ::   ish2
231      REAL(wp),DIMENSION(jpj*jpk) ::   zwork        ! temporary vector for mpp_sum
232      REAL(wp),DIMENSION(jpj,jpk) ::   p_fval       ! return function value
233      !!--------------------------------------------------------------------
234      !
235      p_fval(:,:) = 0.e0
236      IF (PRESENT (bmask)) THEN
237         DO jk = 1, jpkm1
238            DO jj = 2, jpjm1
239               DO ji = fs_2, fs_jpim1   ! Vector opt.
240                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk)                  &
241                     &                          * e1t(ji,jj) * fse3t(ji,jj,jk)   &
242                     &                          * tmask_i(ji,jj)                 &
243                     &                          * bmask(ji,jj)
244               END DO
245            END DO
246         END DO
247      ELSE
248         DO jk = 1, jpkm1
249            DO jj = 2, jpjm1
250               DO ji = fs_2, fs_jpim1   ! Vector opt.
251                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk)                  &
252                     &                          * e1t(ji,jj) * fse3t(ji,jj,jk)   &
253                     &                          * tmask_i(ji,jj)
254               END DO
255            END DO
256         END DO
257      END IF
258      p_fval(:,:) = p_fval(:,:) * 0.5
259      IF(lk_mpp) THEN
260         ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk
261         zwork(:)= RESHAPE( p_fval, ish )
262         CALL mpp_sum( zwork, jpj*jpk )
263         p_fval(:,:)= RESHAPE(zwork,ish2)
264      END IF
265      !
266   END FUNCTION ptr_tjk
267
268
269   SUBROUTINE dia_ptr( kt )
270      !!----------------------------------------------------------------------
271      !!                  ***  ROUTINE dia_ptr  ***
272      !!----------------------------------------------------------------------
273      INTEGER, INTENT(in) ::   kt   ! ocean time step index
274      !!
275      INTEGER  ::   jk, jj, ji               ! dummy loop
276      REAL(wp) ::   zsverdrup,  &              ! conversion from m3/s to Sverdrup
277         &          zpwatt,     &              ! conversion from W    to PW
278         &          zggram                     ! conversion from g    to Pg
279      REAL(wp), DIMENSION(jpi,jpj,jpk) :: vt, vs
280      !!----------------------------------------------------------------------
281
282      IF( kt == nit000 .OR. MOD( kt - nit000 + 1, nf_ptr ) == 0 )   THEN
283
284         zsverdrup = 1.e-6
285         zpwatt    = 1.e-15
286         zggram    = 1.e-6
287         
288         IF ( ln_diaznl ) THEN
289            ! "zonal" mean temperature and salinity at V-points
290            tn_jk_glo(:,:) = ptr_tjk( tn(:,:,:) ) * surf_jk_r_glo(:,:)
291            sn_jk_glo(:,:) = ptr_tjk( sn(:,:,:) ) * surf_jk_r_glo(:,:)
292           
293            IF (ln_subbas) THEN
294               tn_jk_atl(:,:) = ptr_tjk( tn(:,:,:), abasin(:,:) ) * surf_jk_r_atl(:,:)
295               sn_jk_atl(:,:) = ptr_tjk( sn(:,:,:), abasin(:,:) ) * surf_jk_r_atl(:,:)
296               tn_jk_pac(:,:) = ptr_tjk( tn(:,:,:), pbasin(:,:) ) * surf_jk_r_pac(:,:)
297               sn_jk_pac(:,:) = ptr_tjk( sn(:,:,:), pbasin(:,:) ) * surf_jk_r_pac(:,:)
298               tn_jk_ind(:,:) = ptr_tjk( tn(:,:,:), ibasin(:,:) ) * surf_jk_r_ind(:,:)
299               sn_jk_ind(:,:) = ptr_tjk( sn(:,:,:), ibasin(:,:) ) * surf_jk_r_ind(:,:)
300               tn_jk_ipc(:,:) = ptr_tjk( tn(:,:,:), dbasin(:,:) ) * surf_jk_r_ipc(:,:)
301               sn_jk_ipc(:,:) = ptr_tjk( sn(:,:,:), dbasin(:,:) ) * surf_jk_r_ipc(:,:)
302            ENDIF
303         ENDIF
304         
305         !--------------------------------------------------------
306         ! overturning calculation:
307         
308         ! horizontal integral and vertical dz
309
310#if defined key_traldf_eiv
311         v_msf_glo(:,:) = ptr_vjk( vn(:,:,:)+v_eiv(:,:,:) ) 
312         IF( ln_subbas .AND. ln_diaznl ) THEN
313            v_msf_atl(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), abasin(:,:)*sbasin(:,:) ) 
314            v_msf_pac(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), pbasin(:,:)*sbasin(:,:) ) 
315            v_msf_ind(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), ibasin(:,:)*sbasin(:,:) ) 
316            v_msf_ipc(:,:) = ptr_vjk( vn (:,:,:)+v_eiv(:,:,:), dbasin(:,:)*sbasin(:,:) ) 
317         ENDIF
318#else
319         v_msf_glo(:,:) = ptr_vjk( vn(:,:,:) ) 
320         IF( ln_subbas .AND. ln_diaznl ) THEN
321            v_msf_atl(:,:) = ptr_vjk( vn (:,:,:), abasin(:,:)*sbasin(:,:) ) 
322            v_msf_pac(:,:) = ptr_vjk( vn (:,:,:), pbasin(:,:)*sbasin(:,:) ) 
323            v_msf_ind(:,:) = ptr_vjk( vn (:,:,:), ibasin(:,:)*sbasin(:,:) ) 
324            v_msf_ipc(:,:) = ptr_vjk( vn (:,:,:), dbasin(:,:)*sbasin(:,:) ) 
325         ENDIF
326#endif
327         
328#if defined key_diaeiv
329         v_msf_eiv(:,:) = ptr_vjk( v_eiv(:,:,:) ) 
330#endif
331         
332         ! "Meridional" Stream-Function
333         DO jk = 2,jpk 
334            v_msf_glo(:,jk) = v_msf_glo(:,jk-1) + v_msf_glo(:,jk)
335         END DO
336         v_msf_glo(:,:) = v_msf_glo(:,:) * zsverdrup
337#if defined key_diaeiv
338         ! Bolus "Meridional" Stream-Function
339         DO jk = 2,jpk
340            v_msf_eiv(:,jk) = v_msf_eiv(:,jk-1) + v_msf_eiv(:,jk)
341         END DO
342         v_msf_eiv(:,:) = v_msf_eiv(:,:) * zsverdrup
343#endif
344         !
345         IF( ln_subbas .AND. ln_diaznl ) THEN
346            DO jk = 2,jpk 
347               v_msf_atl(:,jk) = v_msf_atl(:,jk-1) + v_msf_atl(:,jk)
348               v_msf_pac(:,jk) = v_msf_pac(:,jk-1) + v_msf_pac(:,jk)
349               v_msf_ind(:,jk) = v_msf_ind(:,jk-1) + v_msf_ind(:,jk)
350               v_msf_ipc(:,jk) = v_msf_ipc(:,jk-1) + v_msf_ipc(:,jk)
351            END DO
352            v_msf_atl(:,:) = v_msf_atl(:,:) * zsverdrup
353            v_msf_pac(:,:) = v_msf_pac(:,:) * zsverdrup
354            v_msf_ind(:,:) = v_msf_ind(:,:) * zsverdrup
355            v_msf_ipc(:,:) = v_msf_ipc(:,:) * zsverdrup
356         ENDIF
357           
358         ! Transports
359         ! T times V on T points (include bolus velocities)
360#if defined key_traldf_eiv 
361         DO jj = 1, jpj
362            DO ji = 1, jpi
363               vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + u_eiv(ji,jj,:) + u_eiv(ji,jj-1,:) )*0.5
364               vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) + v_eiv(ji,jj,:) + v_eiv(ji,jj-1,:) )*0.5
365            END DO
366         END DO
367#else
368         DO jj = 1, jpj
369            DO ji = 1, jpi
370               vt(ji,jj,:) = tn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:) )*0.5
371               vs(ji,jj,:) = sn(ji,jj,:) * ( vn(ji,jj,:) + vn(ji,jj-1,:)  )*0.5
372            END DO
373         END DO
374#endif
375         
376         ht_glo(:) = SUM( ptr_vjk( vt(:,:,:)), 2 )
377         st_glo(:) = SUM( ptr_vjk( vs(:,:,:)), 2 )
378         
379         IF ( ln_subbas ) THEN
380            ht_atl(:) = SUM( ptr_vjk( vt (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
381            ht_pac(:) = SUM( ptr_vjk( vt (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
382            ht_ind(:) = SUM( ptr_vjk( vt (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
383            ht_ipc(:) = SUM( ptr_vjk( vt (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
384            st_atl(:) = SUM( ptr_vjk( vs (:,:,:), abasin(:,:)*sbasin(:,:)), 2 )
385            st_pac(:) = SUM( ptr_vjk( vs (:,:,:), pbasin(:,:)*sbasin(:,:)), 2 )
386            st_ind(:) = SUM( ptr_vjk( vs (:,:,:), ibasin(:,:)*sbasin(:,:)), 2 )
387            st_ipc(:) = SUM( ptr_vjk( vs (:,:,:), dbasin(:,:)*sbasin(:,:)), 2 )
388         ENDIF
389         
390         ! poleward tracer transports:
391         ! overturning components:
392         pht_ove(:) = SUM( v_msf_glo(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
393         pst_ove(:) = SUM( v_msf_glo(:,:) * sn_jk_glo(:,:), 2 )   ! SUM over jk
394#if defined key_diaeiv
395         pht_eiv(:) = SUM( v_msf_eiv(:,:) * tn_jk_glo(:,:), 2 )   ! SUM over jk
396         pst_eiv(:) = SUM( v_msf_eiv(:,:) * sn_jk_glo(:,:), 2 )   ! SUM over jk
397#endif
398
399         ! conversion in PW and G g
400         zpwatt = zpwatt * rau0 * rcp
401         pht_adv(:) = pht_adv(:) * zpwatt 
402         pht_ove(:) = pht_ove(:) * zpwatt
403         pht_ldf(:) = pht_ldf(:) * zpwatt
404         pst_adv(:) = pst_adv(:) * zggram
405         pst_ove(:) = pst_ove(:) * zggram
406         pst_ldf(:) = pst_ldf(:) * zggram
407#if defined key_diaeiv
408         pht_eiv(:) = pht_eiv(:) * zpwatt
409         pst_eiv(:) = pst_eiv(:) * zggram
410#endif
411         IF( ln_subbas ) THEN
412            ht_atl(:) = ht_atl(:) * zpwatt
413            ht_pac(:) = ht_pac(:) * zpwatt
414            ht_ind(:) = ht_ind(:) * zpwatt
415            ht_ipc(:) = ht_ipc(:) * zpwatt
416            st_atl(:) = st_atl(:) * zggram 
417            st_pac(:) = st_pac(:) * zggram
418            st_ind(:) = st_ind(:) * zggram
419            st_ipc(:) = st_ipc(:) * zggram
420         ENDIF
421
422         ! outputs
423         CALL dia_ptr_wri( kt )
424
425      ENDIF
426
427      ! Close the file
428      IF( kt == nitend ) CALL histclo( numptr )
429      !
430   END SUBROUTINE dia_ptr
431
432
433   SUBROUTINE dia_ptr_init
434      !!----------------------------------------------------------------------
435      !!                  ***  ROUTINE dia_ptr_init  ***
436      !!                   
437      !! ** Purpose :   Initialization, namelist read
438      !!----------------------------------------------------------------------
439      NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, nf_ptr, nf_ptr_wri
440      INTEGER ::  inum       ! temporary logical unit
441      !!----------------------------------------------------------------------
442
443      ! Read Namelist namptr : poleward transport parameters
444      REWIND ( numnam )
445      READ   ( numnam, namptr )
446
447      ! Control print
448      IF(lwp) THEN
449         WRITE(numout,*)
450         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
451         WRITE(numout,*) '~~~~~~~~~~~~'
452         WRITE(numout,*) '          Namelist namptr : set ptr parameters'
453         WRITE(numout,*) '             Switch for ptr diagnostic (T) or not (F) ln_diaptr = ', ln_diaptr
454         WRITE(numout,*) '             Atla/Paci/Ind basins computation         ln_subbas = ', ln_subbas
455         WRITE(numout,*) '             Frequency of computation                    nf_ptr = ', nf_ptr
456         WRITE(numout,*) '             Frequency of outputs                    nf_ptr_wri = ', nf_ptr_wri
457      ENDIF
458
459      IF( ln_subbas ) THEN                ! load sub-basin mask
460         CALL iom_open( 'subbasins', inum )
461         CALL iom_get( inum, jpdom_data, 'atlmsk', abasin )      ! Atlantic basin
462         CALL iom_get( inum, jpdom_data, 'pacmsk', pbasin )      ! Pacific basin
463         CALL iom_get( inum, jpdom_data, 'indmsk', ibasin )      ! Indian basin
464         CALL iom_close( inum )
465         dbasin(:,:) = MAX ( pbasin(:,:), ibasin(:,:) )
466         sbasin(:,:) = tmask (:,:,1)
467         WHERE ( gphit (:,:) < -30.e0) sbasin(:,:) = 0.e0
468         IF (lwp) THEN
469            WRITE (numout,*) 'Global ', COUNT (tmask(:,:,1) .GT. 0.5e0 )
470            WRITE (numout,*) 'Atl ', COUNT ( abasin (:,:) .GT. 0.5e0 )
471            WRITE (numout,*) 'Pac ', COUNT ( pbasin (:,:) .GT. 0.5e0 )
472            WRITE (numout,*) 'Ind ', COUNT ( ibasin (:,:) .GT. 0.5e0 )
473            WRITE (numout,*) 'Ipc ', COUNT ( dbasin (:,:) .GT. 0.5e0 )
474            WRITE (numout,*) '30S ', COUNT ( sbasin (:,:) .GT. 0.5e0 )
475         END IF
476      ENDIF
477     
478      ! inverse of the ocean "zonal" v-point section
479      surf_jk_glo(:,:) = ptr_tjk( tmask(:,:,:) )
480      surf_jk_r_glo(:,:) = 0.e0
481      WHERE( surf_jk_glo(:,:) /= 0.e0 )   surf_jk_r_glo(:,:) = 1.e0 / surf_jk_glo(:,:)
482     
483      IF (ln_subbas) THEN
484         surf_jk_atl(:,:) = ptr_tjk( tmask (:,:,:), abasin(:,:) )
485         surf_jk_r_atl(:,:) = 0.e0
486         WHERE( surf_jk_atl(:,:) /= 0.e0 )   surf_jk_r_atl(:,:) = 1.e0 / surf_jk_atl(:,:)
487         !
488         surf_jk_pac(:,:) = ptr_tjk( tmask (:,:,:), pbasin(:,:) )
489         surf_jk_r_pac(:,:) = 0.e0
490         WHERE( surf_jk_pac(:,:) /= 0.e0 )   surf_jk_r_pac(:,:) = 1.e0 / surf_jk_pac(:,:)
491         !
492         surf_jk_ind(:,:) = ptr_tjk( tmask (:,:,:), ibasin(:,:) )
493         surf_jk_r_ind(:,:) = 0.e0
494         WHERE( surf_jk_ind(:,:) /= 0.e0 )   surf_jk_r_ind(:,:) = 1.e0 / surf_jk_ind(:,:)
495         !
496         surf_jk_ipc(:,:) = ptr_tjk( tmask (:,:,:), dbasin(:,:) )
497         surf_jk_r_ipc(:,:) = 0.e0
498         WHERE( surf_jk_ipc(:,:) /= 0.e0 )   surf_jk_r_ipc(:,:) = 1.e0 / surf_jk_ipc(:,:)
499      END IF
500     
501   END SUBROUTINE dia_ptr_init
502
503
504   SUBROUTINE dia_ptr_wri( kt )
505      !!---------------------------------------------------------------------
506      !!                ***  ROUTINE dia_ptr_wri  ***
507      !!
508      !! ** Purpose :   output of poleward fluxes
509      !!
510      !! ** Method  :   NetCDF file
511      !!----------------------------------------------------------------------
512      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
513      !!
514      INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw
515      INTEGER, SAVE :: ndim  , ndim_atl     , ndim_pac     , ndim_ind     , ndim_ipc
516      INTEGER, SAVE ::         ndim_atl_30  , ndim_pac_30  , ndim_ind_30  , ndim_ipc_30
517      INTEGER, SAVE :: ndim_h, ndim_h_atl_30, ndim_h_pac_30, ndim_h_ind_30, ndim_h_ipc_30
518      INTEGER, SAVE, DIMENSION (jpj*jpk) :: ndex  , ndex_atl     , ndex_pac     , ndex_ind     , ndex_ipc
519      INTEGER, SAVE, DIMENSION (jpj*jpk) ::         ndex_atl_30  , ndex_pac_30  , ndex_ind_30  , ndex_ipc_30
520      INTEGER, SAVE, DIMENSION (jpj)     :: ndex_h, ndex_h_atl_30, ndex_h_pac_30, ndex_h_ind_30, ndex_h_ipc_30
521
522      CHARACTER (len=40)       ::   clhstnam, clop, clop_once         ! temporary names
523      INTEGER                  ::   iline, it, ji, jj, jk             !
524      REAL(wp)                 ::   zsto, zout, zdt, zjulian   ! temporary scalars
525      REAL(wp), DIMENSION(jpj) ::   zphi, zfoo
526      REAL(wp), DIMENSION(jpj,jpk) :: z_1
527      !!----------------------------------------------------------------------
528
529      ! define time axis
530      it = kt - nit000 + 1
531
532      ! Initialization
533      ! --------------
534      IF( kt == nit000 ) THEN
535
536         zdt = rdt
537         IF( nacc == 1 ) zdt = rdtmin
538
539         ! Reference latitude
540         ! ------------------
541         !                                           ! =======================
542         IF( cp_cfg == "orca" ) THEN                 !   ORCA configurations
543            !                                        ! =======================
544
545            IF( jp_cfg == 05  )   iline = 192   ! i-line that passes near the North Pole
546            IF( jp_cfg == 025 )   iline = 384   ! i-line that passes near the North Pole
547            IF( jp_cfg == 2   )   iline =  48   ! i-line that passes near the North Pole
548            IF( jp_cfg == 4   )   iline =  24   ! i-line that passes near the North Pole
549            zphi(:) = 0.e0
550            DO ji = mi0(iline), mi1(iline) 
551               zphi(:) = gphiv(ji,:)         ! if iline is in the local domain
552               ! correct highest latitude for some configurations - will work if domain is parallelized in J ?
553               IF( jp_cfg == 05 ) THEN
554                  DO jj = mj0(jpj), mj1(jpj) 
555                     zphi(jj) = zphi(jpjm1) + (zphi(jpjm1)-zphi(jpj-2))/2.
556                     zphi(jj) = MIN( zphi(jj), 90.)
557                  END DO
558               END IF
559               IF( jp_cfg == 2 ) THEN
560                  DO jj = mj0(jpj-1), mj1(jpj-1) 
561                     zphi(jj  ) = 88.5e0
562                  END DO
563                  DO jj = mj0(jpj  ), mj1(jpj  ) 
564                     zphi(jj  ) = 89.5e0
565                  END DO
566               END IF
567            END DO
568            ! provide the correct zphi to all local domains
569            IF( lk_mpp )   CALL mpp_sum( zphi, jpj )       
570
571            !                                        ! =======================
572         ELSE                                        !   OTHER configurations
573            !                                        ! =======================
574            zphi(:) = gphiv(1,:)             ! assume lat/lon coordinate, select the first i-line
575            !
576         ENDIF
577
578         ! OPEN netcdf file
579         ! ----------------
580         ! Define frequency of output and means
581         zsto = nf_ptr * zdt
582         IF( ln_mskland )   THEN    ! put 1.e+20 on land (very expensive!!)
583            clop      = "ave(only(x))"
584            clop_once = "once(only(x))"
585         ELSE                       ! no use of the mask value (require less cpu time)
586            clop      = "ave(x)"       
587            clop_once = "once"
588         ENDIF
589
590         zout = nf_ptr_wri * zdt
591         zfoo(:) = 0.e0
592
593         ! Compute julian date from starting date of the run
594
595         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
596         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
597
598         CALL dia_nam( clhstnam, nf_ptr_wri, 'diaptr' )
599         IF(lwp)WRITE( numout,*)" Name of diaptr NETCDF file ",clhstnam
600
601         ! Horizontal grid : zphi()
602         CALL histbeg(clhstnam, 1, zfoo, jpj, zphi,   &
603            1, 1, 1, jpj, 0, zjulian, zdt, nhoridz, numptr, domain_id=nidom )
604         ! Vertical grids : gdept_0, gdepw_0
605         CALL histvert( numptr, "deptht", "Vertical T levels",   &
606            "m", jpk, gdept_0, ndepidzt )
607         CALL histvert( numptr, "depthw", "Vertical W levels",   &
608            "m", jpk, gdepw_0, ndepidzw )
609
610         !
611         CALL wheneq ( jpj*jpk, MIN(surf_jk_glo(:,:), 1.e0), 1, 1., ndex  , ndim  )      ! Lat-Depth
612         CALL wheneq ( jpj    , MIN(surf_jk_glo(:,1), 1.e0), 1, 1., ndex_h, ndim_h )     ! Lat
613         WRITE (numout,*) 'ndex : ', COUNT(surf_jk_glo>0.), COUNT(surf_jk_r_glo>0.), SUM(ndex), SUM(ndex_h), ndim, ndim_h 
614         WRITE (numout,*) 'ndex : ', SUM(ndex), SUM(ndex_h), ndim, ndim_h 
615
616         IF (ln_subbas) THEN
617            z_1 (:,1) = 1.0e0
618            WHERE ( gphit (jpi/2,:) .LT. -30 ) z_1 (:,1) = 0.e0
619            DO jk = 2, jpk
620               z_1 (:,jk) = z_1 (:,1)
621            END DO
622
623            CALL wheneq ( jpj*jpk, MIN(surf_jk_atl(:,:)         , 1.e0), 1, 1., ndex_atl     , ndim_atl      ) ! Lat-Depth
624            CALL wheneq ( jpj*jpk, MIN(surf_jk_atl(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_atl_30  , ndim_atl_30   ) ! Lat-Depth
625            CALL wheneq ( jpj    , MIN(surf_jk_atl(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_atl_30, ndim_h_atl_30 ) ! Lat
626            WRITE (numout,*) 'ndex atl : ', SUM(ndex_atl), SUM(ndex_atl_30), SUM(ndex_h_atl_30),&
627               & ndim_atl, ndim_atl_30, ndim_h_atl_30
628
629            CALL wheneq ( jpj*jpk, MIN(surf_jk_pac(:,:)         , 1.e0), 1, 1., ndex_pac     , ndim_pac      ) ! Lat-Depth
630            CALL wheneq ( jpj*jpk, MIN(surf_jk_pac(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_pac_30  , ndim_pac_30   ) ! Lat-Depth
631            CALL wheneq ( jpj    , MIN(surf_jk_pac(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_pac_30, ndim_h_pac_30 ) ! Lat
632            WRITE (numout,*) 'ndex pac : ', SUM(ndex_pac), SUM(ndex_pac_30), SUM(ndex_h_pac_30),&
633               & ndim_pac, ndim_pac_30, ndim_h_pac_30
634
635            CALL wheneq ( jpj*jpk, MIN(surf_jk_ind(:,:)         , 1.e0), 1, 1., ndex_ind     , ndim_ind      ) ! Lat-Depth
636            CALL wheneq ( jpj*jpk, MIN(surf_jk_ind(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_ind_30  , ndim_ind_30   ) ! Lat-Depth
637            CALL wheneq ( jpj    , MIN(surf_jk_ind(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_ind_30, ndim_h_ind_30 ) ! Lat
638            WRITE (numout,*) 'ndex ind : ', SUM(ndex_ind), SUM(ndex_ind_30), SUM(ndex_h_ind_30),&
639               & ndim_ind, ndim_ind_30, ndim_h_ind_30
640
641            CALL wheneq ( jpj*jpk, MIN(surf_jk_ipc(:,:)         , 1.e0), 1, 1., ndex_ipc     , ndim_ipc      ) ! Lat-Depth
642            CALL wheneq ( jpj*jpk, MIN(surf_jk_ipc(:,:)*z_1(:,:), 1.e0), 1, 1., ndex_ipc_30  , ndim_ipc_30   ) ! Lat-Depth
643            CALL wheneq ( jpj    , MIN(surf_jk_ipc(:,1)*z_1(:,1), 1.e0), 1, 1., ndex_h_ipc_30, ndim_h_ipc_30 ) ! Lat
644            WRITE (numout,*) 'ndex ipc : ', SUM(ndex_ipc), SUM(ndex_ipc_30), SUM(ndex_h_ipc_30),&
645               & ndim_ipc, ndim_ipc_30, ndim_h_ipc_30
646
647         ENDIF
648
649         !  Zonal mean T and S
650
651         IF ( ln_diaznl ) THEN
652            CALL histdef( numptr, "zotemglo", "Zonal Mean Temperature","C" ,   &
653               1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
654            CALL histdef( numptr, "zosalglo", "Zonal Mean Salinity","PSU"  ,   &
655               1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
656
657            CALL histdef( numptr, "zosrfglo", "Zonal Mean Surface","m^2"  ,   &
658               1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
659
660            IF (ln_subbas) THEN
661               CALL histdef( numptr, "zotematl", "Zonal Mean Temperature: Atlantic","C" ,   &
662                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
663               CALL histdef( numptr, "zosalatl", "Zonal Mean Salinity: Atlantic","PSU"  ,   &
664                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
665               CALL histdef( numptr, "zosrfatl", "Zonal Mean Surface: Atlantic","m^2"  ,   &
666                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
667
668               CALL histdef( numptr, "zotempac", "Zonal Mean Temperature: Pacific","C" ,   &
669                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
670               CALL histdef( numptr, "zosalpac", "Zonal Mean Salinity: Pacific","PSU"  ,   &
671                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
672               CALL histdef( numptr, "zosrfpac", "Zonal Mean Surface: Pacific","m^2"  ,   &
673                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
674
675               CALL histdef( numptr, "zotemind", "Zonal Mean Temperature: Indian","C" ,   &
676                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
677               CALL histdef( numptr, "zosalind", "Zonal Mean Salinity: Indian","PSU"  ,   &
678                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
679               CALL histdef( numptr, "zosrfind", "Zonal Mean Surface: Indian","m^2"  ,   &
680                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
681
682               CALL histdef( numptr, "zotemipc", "Zonal Mean Temperature: Pacific+Indian","C" ,   &
683                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
684               CALL histdef( numptr, "zosalipc", "Zonal Mean Salinity: Pacific+Indian","PSU"  ,   &
685                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
686               CALL histdef( numptr, "zosrfipc", "Zonal Mean Surface: Pacific+Indian","m^2"  ,   &
687                  1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout )
688            ENDIF
689
690         ENDIF
691
692         !  Meridional Stream-Function (eulerian and bolus)
693
694         CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: Global","Sv" ,   &
695            1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
696         IF( ln_subbas .AND. ln_diaznl ) THEN
697            CALL histdef( numptr, "zomsfatl", "Meridional Stream-Function: Atlantic","Sv" ,   &
698               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
699            CALL histdef( numptr, "zomsfpac", "Meridional Stream-Function: Pacific","Sv" ,   &
700               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
701            CALL histdef( numptr, "zomsfind", "Meridional Stream-Function: Indian","Sv" ,   &
702               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
703            CALL histdef( numptr, "zomsfipc", "Meridional Stream-Function: Indo-Pacific","Sv" ,&
704               1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
705         ENDIF
706
707         !  Heat transport
708
709         CALL histdef( numptr, "sophtadv", "Advective Heat Transport"      ,   &
710            "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
711         CALL histdef( numptr, "sophtldf", "Diffusive Heat Transport"      ,   &
712            "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
713         CALL histdef( numptr, "sophtove", "Overturning Heat Transport"    ,   &
714            "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
715         IF( ln_subbas ) THEN
716            CALL histdef( numptr, "sohtatl", "Heat Transport Atlantic"      ,  &
717               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
718            CALL histdef( numptr, "sohtpac", "Heat Transport Pacific"      ,   &
719               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
720            CALL histdef( numptr, "sohtind", "Heat Transport Indian"      ,     &
721               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
722            CALL histdef( numptr, "sohtipc", "Heat Transport Pacific+Indian"      ,     &
723               "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
724         ENDIF
725
726
727         !  Salt transport
728
729         CALL histdef( numptr, "sopstadv", "Advective Salt Transport"      ,   &
730            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
731         CALL histdef( numptr, "sopstldf", "Diffusive Salt Transport"      ,   &
732            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
733         CALL histdef( numptr, "sopstove", "Overturning Salt Transport"    ,   &
734            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
735
736#if defined key_diaeiv
737         ! Eddy induced velocity
738         CALL histdef( numptr, "zomsfeiv", "Bolus Meridional Stream-Function: global",   &
739            "Sv"      , 1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
740         CALL histdef( numptr, "sophteiv", "Bolus Advective Heat Transport",   &
741            "PW"      , 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
742         CALL histdef( numptr, "sopsteiv", "Bolus Advective Salt Transport",   &
743            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
744#endif
745         IF( ln_subbas ) THEN
746            CALL histdef( numptr, "sostatl", "Salt Transport Atlantic"      ,    &
747               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
748            CALL histdef( numptr, "sostpac", "Salt Transport Pacific"      ,     &
749               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
750            CALL histdef( numptr, "sostind", "Salt Transport Indian"      ,      &
751               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
752            CALL histdef( numptr, "sostipc", "Salt Transport Pacific+Indian",     &
753               "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
754         ENDIF
755
756         CALL histend( numptr )
757
758      ENDIF
759
760      IF( MOD( it, nf_ptr ) == 0 ) THEN
761
762         IF(lwp) THEN
763            WRITE(numout,*)
764            WRITE(numout,*) 'dia_ptr : compute Poleward Transports at time-step : ', kt
765            WRITE(numout,*) '~~~~~~~~'
766            WRITE(numout,*)
767         ENDIF
768
769         IF (ln_diaznl ) THEN
770            CALL histwrite( numptr, "zosrfglo", it, surf_jk_glo , ndim, ndex )
771            CALL histwrite( numptr, "zotemglo", it, tn_jk_glo  , ndim, ndex )
772            CALL histwrite( numptr, "zosalglo", it, sn_jk_glo  , ndim, ndex )
773
774            IF (ln_subbas) THEN
775               CALL histwrite( numptr, "zosrfatl", it, surf_jk_atl, ndim_atl, ndex_atl )
776               CALL histwrite( numptr, "zosrfpac", it, surf_jk_pac, ndim_pac, ndex_pac )
777               CALL histwrite( numptr, "zosrfind", it, surf_jk_ind, ndim_ind, ndex_ind )
778               CALL histwrite( numptr, "zosrfipc", it, surf_jk_ipc, ndim_ipc, ndex_ipc )
779
780               CALL histwrite( numptr, "zotematl", it, tn_jk_atl  , ndim_atl, ndex_atl )
781               CALL histwrite( numptr, "zosalatl", it, sn_jk_atl  , ndim_atl, ndex_atl )
782               CALL histwrite( numptr, "zotempac", it, tn_jk_pac  , ndim_pac, ndex_pac )
783               CALL histwrite( numptr, "zosalpac", it, sn_jk_pac  , ndim_pac, ndex_pac )
784               CALL histwrite( numptr, "zotemind", it, tn_jk_ind  , ndim_ind, ndex_ind )
785               CALL histwrite( numptr, "zosalind", it, sn_jk_ind  , ndim_ind, ndex_ind )
786               CALL histwrite( numptr, "zotemipc", it, tn_jk_ipc  , ndim_ipc, ndex_ipc )
787               CALL histwrite( numptr, "zosalipc", it, sn_jk_ipc  , ndim_ipc, ndex_ipc )
788            END IF
789         ENDIF
790
791         ! overturning outputs:
792         CALL histwrite( numptr, "zomsfglo", it, v_msf_glo, ndim, ndex )
793         IF( ln_subbas .AND. ln_diaznl ) THEN
794            CALL histwrite( numptr, "zomsfatl", it, v_msf_atl , ndim_atl_30, ndex_atl_30 )
795            CALL histwrite( numptr, "zomsfpac", it, v_msf_pac , ndim_pac_30, ndex_pac_30 )
796            CALL histwrite( numptr, "zomsfind", it, v_msf_ind , ndim_ind_30, ndex_ind_30 )
797            CALL histwrite( numptr, "zomsfipc", it, v_msf_ipc , ndim_ipc_30, ndex_ipc_30 )
798         ENDIF
799#if defined key_diaeiv
800         CALL histwrite( numptr, "zomsfeiv", it, v_msf_eiv, ndim  , ndex   )
801#endif
802
803      ENDIF
804
805      ! heat transport outputs:
806      IF( ln_subbas ) THEN
807         CALL histwrite( numptr, "sohtatl", it, ht_atl  , ndim_h_atl_30, ndex_h_atl_30 )
808         CALL histwrite( numptr, "sohtpac", it, ht_pac  , ndim_h_pac_30, ndex_h_pac_30 )
809         CALL histwrite( numptr, "sohtind", it, ht_ind  , ndim_h_ind_30, ndex_h_ind_30 )
810         CALL histwrite( numptr, "sohtipc", it, ht_ipc  , ndim_h_ipc_30, ndex_h_ipc_30 )
811         CALL histwrite( numptr, "sostatl", it, st_atl  , ndim_h_atl_30, ndex_h_atl_30 )
812         CALL histwrite( numptr, "sostpac", it, st_pac  , ndim_h_pac_30, ndex_h_pac_30 )
813         CALL histwrite( numptr, "sostind", it, st_ind  , ndim_h_ind_30, ndex_h_ind_30 )
814         CALL histwrite( numptr, "sostipc", it, st_ipc  , ndim_h_ipc_30, ndex_h_ipc_30 )
815      ENDIF
816
817      CALL histwrite( numptr, "sophtadv", it, pht_adv  , ndim_h, ndex_h )
818      CALL histwrite( numptr, "sophtldf", it, pht_ldf  , ndim_h, ndex_h )
819      CALL histwrite( numptr, "sophtove", it, pht_ove  , ndim_h, ndex_h )
820      CALL histwrite( numptr, "sopstadv", it, pst_adv  , ndim_h, ndex_h )
821      CALL histwrite( numptr, "sopstldf", it, pst_ldf  , ndim_h, ndex_h )
822      CALL histwrite( numptr, "sopstove", it, pst_ove  , ndim_h, ndex_h )
823#if defined key_diaeiv
824      CALL histwrite( numptr, "sophteiv", it, pht_eiv  , ndim_h, ndex_h )
825      CALL histwrite( numptr, "sopsteiv", it, pst_eiv  , ndim_h, ndex_h )
826#endif
827
828      !
829   END SUBROUTINE dia_ptr_wri
830
831   !!======================================================================
832END MODULE diaptr