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

source: branches/2015/dev_r5020_CNRS_DIAPTR/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90 @ 5137

Last change on this file since 5137 was 5137, checked in by cetlod, 9 years ago

dev_r5020_CNRS_DIAPTR: some improvments of poleward heat transport diagnostics before merging with the trunk
  • Property svn:keywords set to Id
File size: 17.9 KB
Line 
1MODULE diaptr
2   !!======================================================================
3   !!                       ***  MODULE  diaptr  ***
4   !! Ocean physics:  Computes meridonal transports and zonal means
5   !!=====================================================================
6   !! History :  1.0  ! 2003-09  (C. Talandier, G. Madec)  Original code
7   !!            2.0  ! 2006-01  (A. Biastoch)  Allow sub-basins computation
8   !!            3.2  ! 2010-03  (O. Marti, S. Flavoni) Add fields
9   !!            3.3  ! 2010-10  (G. Madec)  dynamical allocation
10   !!            3.6  ! 2014-12  (C. Ethe) use of IOM
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dia_ptr      : Poleward Transport Diagnostics module
15   !!   dia_ptr_init : Initialization, namelist read
16   !!   ptr_sjk      : "zonal" mean computation of a field - tracer or flux array
17   !!   ptr_sj       : "zonal" and vertical sum computation of a "meridional" flux array
18   !!                   (Generic interface to ptr_sj_3d, ptr_sj_2d)
19   !!----------------------------------------------------------------------
20   USE oce              ! ocean dynamics and active tracers
21   USE dom_oce          ! ocean space and time domain
22   USE domngb
23   USE phycst           ! physical constants
24   !
25   USE iom              ! IOM library
26   USE in_out_manager   ! I/O manager
27   USE lib_mpp          ! MPP library
28   USE timing           ! preformance summary
29
30   IMPLICIT NONE
31   PRIVATE
32
33   INTERFACE ptr_sj
34      MODULE PROCEDURE ptr_sj_3d, ptr_sj_2d
35   END INTERFACE
36
37   PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines
38   PUBLIC   ptr_sjk        !
39   PUBLIC   dia_ptr_init   ! call in step module
40   PUBLIC   dia_ptr        ! call in step module
41
42   !                                  !!** namelist  namptr  **
43   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf   !: Heat TRansports (adv, diff, overturn.)
44   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf   !: Salt TRansports (adv, diff, overturn.)
45   
46
47   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F)
48   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation
49   INTEGER         ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T)
50   INTEGER         ::   nx, ny 
51
52   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
53   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rau0 x Cp)
54   REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg
55
56   CHARACTER(len=3), ALLOCATABLE, SAVE, DIMENSION(:)     :: clsubb
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
60   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)     :: p_fval1d
61   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: p_fval2d
62
63
64   !! * Substitutions
65#  include "domzgr_substitute.h90"
66#  include "vectopt_loop_substitute.h90"
67   !!----------------------------------------------------------------------
68   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
69   !! $Id$
70   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
71   !!----------------------------------------------------------------------
72CONTAINS
73
74   SUBROUTINE dia_ptr( pvtr )
75      !!----------------------------------------------------------------------
76      !!                  ***  ROUTINE dia_ptr  ***
77      !!----------------------------------------------------------------------
78      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
79      !
80      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
81      REAL(wp) ::   zv, zsfc               ! local scalar
82      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
83      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace
84      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace
85      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace
86      CHARACTER( len = 10 )  :: cl1
87      !!----------------------------------------------------------------------
88      !
89      IF( nn_timing == 1 )   CALL timing_start('dia_ptr')
90
91      !
92      IF( PRESENT( pvtr ) ) THEN
93         IF( iom_use("zomsfglo") ) THEN    ! effective MSF
94            z3d(nx,:,:) = ptr_sjk( pvtr(:,:,:) )  ! zonal cumulative effective transport
95            DO jk = 2, jpkm1 
96              z3d(nx,:,jk) = z3d(nx,:,jk-1) + z3d(nx,:,jk)   ! effective j-Stream-Function (MSF)
97           END DO
98            cl1 = TRIM('zomsf'//clsubb(1) )
99            CALL iom_put( cl1, z3d * rc_sv )
100            DO jn = 2, nptr                                    ! by sub-basins
101               z3d(nx,:,:) =  ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) 
102               DO jk = 2, jpkm1 
103                  z3d(nx,:,jk) = z3d(nx,:,jk-1) + z3d(nx,:,jk)    ! effective j-Stream-Function (MSF)
104               END DO
105               cl1 = TRIM('zomsf'//clsubb(jn) )
106               CALL iom_put( cl1, z3d * rc_sv )
107            END DO
108         ENDIF
109         !
110      ELSE
111         !
112         IF( iom_use("zotemglo") ) THEN    ! i-mean i-k-surface
113            DO jk = 1, jpkm1
114               DO jj = 1, jpj
115                  DO ji = 1, jpi
116                     zsfc = e1t(ji,jj) * fse3t(ji,jj,jk)
117                     zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
118                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc
119                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc
120                  ENDDO
121               ENDDO
122            ENDDO
123            DO jn = 1, nptr
124               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
125               cl1 = TRIM('zosrf'//clsubb(jn) )
126               CALL iom_put( cl1, zmask )
127               !
128               z3d(nx,:,:) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) &
129                  &            / MAX( zmask(1,:,:), 10.e-15 )
130               cl1 = TRIM('zotem'//clsubb(jn) )
131               CALL iom_put( cl1, z3d )
132               !
133               z3d(nx,:,:) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) &
134                  &            / MAX( zmask(1,:,:), 10.e-15 )
135               cl1 = TRIM('zosal'//clsubb(jn) )
136               CALL iom_put( cl1, z3d )
137            END DO
138         ENDIF
139         !
140         !                                ! Advective and diffusive heat and salt transport
141         IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN   
142            z2d(nx,:) = htr_adv(:) * rc_pwatt        !  (conversion in PW)
143            cl1 = 'sophtadv'                 
144            !IF(lwp) write(numout,*) ' sophtadv : ',(z2d(nx,jj), jj=1,jpj)
145            CALL iom_put( TRIM(cl1), z2d )
146            z2d(nx,:) = str_adv(:) * rc_ggram        ! (conversion in Gg)
147            cl1 = 'sopstadv'
148            CALL iom_put( TRIM(cl1), z2d )
149         ENDIF
150         !
151         IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN   
152            z2d(nx,:) = htr_ldf(:) * rc_pwatt        !  (conversion in PW)
153            cl1 = 'sophtldf'
154            CALL iom_put( TRIM(cl1), z2d )
155            z2d(nx,:) = str_ldf(:) * rc_ggram        !  (conversion in Gg)
156            cl1 = 'sopstldf'
157            CALL iom_put( TRIM(cl1), z2d )
158         ENDIF
159         !
160      ENDIF
161      !
162      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr')
163      !
164   END SUBROUTINE dia_ptr
165
166
167   SUBROUTINE dia_ptr_init
168      !!----------------------------------------------------------------------
169      !!                  ***  ROUTINE dia_ptr_init  ***
170      !!                   
171      !! ** Purpose :   Initialization, namelist read
172      !!----------------------------------------------------------------------
173      INTEGER ::  jn           ! local integers
174      INTEGER ::  inum, ierr   ! local integers
175      INTEGER ::  ios          ! Local integer output status for namelist read
176      !!
177      NAMELIST/namptr/ ln_diaptr, ln_subbas
178      !!----------------------------------------------------------------------
179
180      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
181      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901)
182901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist', lwp )
183
184      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport
185      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 )
186902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp )
187      IF(lwm) WRITE ( numond, namptr )
188
189      IF(lwp) THEN                     ! Control print
190         WRITE(numout,*)
191         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
192         WRITE(numout,*) '~~~~~~~~~~~~'
193         WRITE(numout,*) '   Namelist namptr : set ptr parameters'
194         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr
195         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas
196      ENDIF
197
198      IF( ln_diaptr ) THEN 
199         !
200         CALL dom_ngb( 180., 90., nx, ny, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots)
201         !
202         IF( ln_subbas ) THEN
203            nptr = 5            ! Global, Atlantic, Pacific, Indian, Indo-Pacific
204            ALLOCATE( clsubb(nptr) )
205            clsubb(1) = 'glo' ;  clsubb(2) = 'atl'  ;  clsubb(3) = 'pac'  ;  clsubb(4) = 'ind'  ;  clsubb(5) = 'ipc'
206         ELSE               
207            nptr = 1       ! Global only
208            ALLOCATE( clsubb(nptr) )
209            clsubb(1) = 'glo' 
210         ENDIF
211
212         !                                      ! allocate dia_ptr arrays
213         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
214
215         rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt
216
217         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
218
219         IF( ln_subbas ) THEN                ! load sub-basin mask
220            CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  )
221            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
222            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
223            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
224            CALL iom_close( inum )
225            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
226            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
227            ELSE WHERE                     ;   btm30(:,:) = ssmask(:,:)
228            END WHERE
229         ENDIF
230   
231         btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
232     
233         DO jn = 1, nptr
234            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
235         END DO
236
237         ! Initialise arrays to zero because diatpr is called before they are first calculated
238         ! Note that this means diagnostics will not be exactly correct when model run is restarted.
239         htr_adv(:) = 0._wp  ;  str_adv(:) =  0._wp 
240         htr_ldf(:) = 0._wp  ;  str_ldf(:) =  0._wp 
241         !
242      ENDIF 
243      !
244   END SUBROUTINE dia_ptr_init
245
246
247   FUNCTION dia_ptr_alloc()
248      !!----------------------------------------------------------------------
249      !!                    ***  ROUTINE dia_ptr_alloc  ***
250      !!----------------------------------------------------------------------
251      INTEGER               ::   dia_ptr_alloc   ! return value
252      INTEGER, DIMENSION(3) ::   ierr
253      !!----------------------------------------------------------------------
254      ierr(:) = 0
255      !
256      ALLOCATE( btmsk(jpi,jpj,nptr) ,           &
257         &      htr_adv(jpj) , str_adv(jpj) ,   &
258         &      htr_ldf(jpj) , str_ldf(jpj) , STAT=ierr(1)  )
259         !
260      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2))
261      !
262      ALLOCATE( btm30(jpi,jpj), STAT=ierr(3)  )
263
264         !
265      dia_ptr_alloc = MAXVAL( ierr )
266      IF(lk_mpp)   CALL mpp_sum( dia_ptr_alloc )
267      !
268   END FUNCTION dia_ptr_alloc
269
270
271   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval )
272      !!----------------------------------------------------------------------
273      !!                    ***  ROUTINE ptr_sj_3d  ***
274      !!
275      !! ** Purpose :   i-k sum computation of a j-flux array
276      !!
277      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
278      !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
279      !!
280      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
281      !!----------------------------------------------------------------------
282      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pva   ! mask flux array at V-point
283      REAL(wp), INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
284      !
285      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
286      INTEGER                  ::   ijpj         ! ???
287      REAL(wp), POINTER, DIMENSION(:) :: p_fval  ! function value
288      !!--------------------------------------------------------------------
289      !
290      p_fval => p_fval1d
291
292      ijpj = jpj
293      p_fval(:) = 0._wp
294      IF( PRESENT( pmsk ) ) THEN
295         DO jk = 1, jpkm1
296            DO jj = 2, jpjm1
297               DO ji = fs_2, fs_jpim1   ! Vector opt.
298                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) * pmsk(ji,jj)
299               END DO
300            END DO
301         END DO
302      ELSE
303         DO jk = 1, jpkm1
304            DO jj = 2, jpjm1
305               DO ji = fs_2, fs_jpim1   ! Vector opt.
306                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) 
307               END DO
308            END DO
309         END DO
310      ENDIF
311#if defined key_mpp_mpi
312      IF(lk_mpp)   CALL mpp_sum( p_fval, ijpj, ncomm_znl)
313#endif
314      !
315   END FUNCTION ptr_sj_3d
316
317
318   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval )
319      !!----------------------------------------------------------------------
320      !!                    ***  ROUTINE ptr_sj_2d  ***
321      !!
322      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array
323      !!
324      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
325      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
326      !!
327      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
328      !!----------------------------------------------------------------------
329      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pva   ! mask flux array at V-point
330      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
331      !
332      INTEGER                  ::   ji,jj       ! dummy loop arguments
333      INTEGER                  ::   ijpj        ! ???
334      REAL(wp), POINTER, DIMENSION(:) :: p_fval ! function value
335      !!--------------------------------------------------------------------
336      !
337      p_fval => p_fval1d
338
339      ijpj = jpj
340      p_fval(:) = 0._wp
341      IF( PRESENT( pmsk ) ) THEN
342         DO jj = 2, jpjm1
343            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
344               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) * pmsk(ji,jj)
345            END DO
346         END DO
347      ELSE
348         DO jj = 2, jpjm1
349            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
350               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj)
351            END DO
352         END DO
353      ENDIF
354#if defined key_mpp_mpi
355      CALL mpp_sum( p_fval, ijpj, ncomm_znl )
356#endif
357      !
358   END FUNCTION ptr_sj_2d
359
360
361   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
362      !!----------------------------------------------------------------------
363      !!                    ***  ROUTINE ptr_sjk  ***
364      !!
365      !! ** Purpose :   i-sum computation of an array
366      !!
367      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
368      !!
369      !! ** Action  : - p_fval: i-mean poleward flux of pva
370      !!----------------------------------------------------------------------
371      !!
372      IMPLICIT none
373      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pta    ! mask flux array at V-point
374      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask
375      !!
376      INTEGER                           :: ji, jj, jk ! dummy loop arguments
377      REAL(wp), POINTER, DIMENSION(:,:) :: p_fval     ! return function value
378#if defined key_mpp_mpi
379      INTEGER, DIMENSION(1) ::   ish
380      INTEGER, DIMENSION(2) ::   ish2
381      INTEGER               ::   ijpjjpk
382      REAL(wp), DIMENSION(jpj*jpk) ::   zwork    ! mask flux array at V-point
383#endif
384      !!--------------------------------------------------------------------
385      !
386      p_fval => p_fval2d
387
388      p_fval(:,:) = 0._wp
389      !
390      IF( PRESENT( pmsk ) ) THEN
391         DO jk = 1, jpkm1
392            DO jj = 2, jpjm1
393!!gm here, use of tmask_i  ==> no need of loop over nldi, nlei....
394               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
395                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj)
396               END DO
397            END DO
398         END DO
399      ELSE
400         DO jk = 1, jpkm1
401            DO jj = 2, jpjm1
402               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
403                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * tmask_i(ji,jj)
404               END DO
405            END DO
406         END DO
407      END IF
408      !
409#if defined key_mpp_mpi
410      ijpjjpk = jpj*jpk
411      ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk
412      zwork(1:ijpjjpk) = RESHAPE( p_fval, ish )
413      CALL mpp_sum( zwork, ijpjjpk, ncomm_znl )
414      p_fval(:,:) = RESHAPE( zwork, ish2 )
415#endif
416      !
417   END FUNCTION ptr_sjk
418
419
420   !!======================================================================
421END MODULE diaptr
Note: See TracBrowser for help on using the repository browser.