source: CONFIG/UNIFORM/v6/IPSLCM6CHT/SOURCES/NEMO/diaptr.F90 @ 4460

Last change on this file since 4460 was 2456, checked in by acosce, 9 years ago

Add new configuration IPSLCM6CHT

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