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 @ 5023

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

branch 2015/dev_r5020_CNRS_DIAPTR : Poleward TRansports diagnostics using XIOS, see ticket #1435
  • Property svn:keywords set to Id
File size: 17.5 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_alloc   ! call in opa 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 ::   l_diaptr = .TRUE.    !  Poleward transport flag (T) or not (F)
48   LOGICAL         ::   l_subbas = .FALSE.   !  Atlantic/Pacific/Indian basins calculation
49   LOGICAL         ::   linit    = .TRUE.    !  initialization flag (set to false after the 1st call)
50   INTEGER         ::   nptr                 ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T)
51   INTEGER         ::   nx, ny 
52
53   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
54   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rau0 x Cp)
55   REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg
56
57   CHARACTER(len=3), ALLOCATABLE, SAVE, DIMENSION(:)     :: clsubb
58   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
59   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   :: btm30   ! mask out Southern Ocean (=0 south of 30°S)
60
61   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)     :: p_fval1d
62   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: p_fval2d
63
64
65   !! * Substitutions
66#  include "domzgr_substitute.h90"
67#  include "vectopt_loop_substitute.h90"
68   !!----------------------------------------------------------------------
69   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
70   !! $Id$
71   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
72   !!----------------------------------------------------------------------
73CONTAINS
74
75   SUBROUTINE dia_ptr( pvtr )
76      !!----------------------------------------------------------------------
77      !!                  ***  ROUTINE dia_ptr  ***
78      !!----------------------------------------------------------------------
79      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
80      !
81      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
82      REAL(wp) ::   zv, zsfc               ! local scalar
83      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
84      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace
85      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace
86      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace
87      CHARACTER( len = 10 )  :: cl1
88      !!----------------------------------------------------------------------
89      !
90      IF( nn_timing == 1 )   CALL timing_start('dia_ptr')
91
92      IF( linit )  THEN
93         CALL dia_ptr_init
94         linit = .FALSE.
95      ENDIF
96      !
97      IF( PRESENT( pvtr ) ) THEN
98         IF( iom_use("zomsfglo") ) THEN    ! effective MSF
99            z3d(nx,:,:) = ptr_sjk( pvtr(:,:,:) )  ! zonal cumulative effective transport
100            DO jk = 2, jpkm1 
101              z3d(nx,:,jk) = z3d(nx,:,jk-1) + z3d(nx,:,jk)   ! effective j-Stream-Function (MSF)
102           END DO
103            cl1 = TRIM('zomsf'//clsubb(1) )
104            CALL iom_put( cl1, z3d * rc_sv )
105            DO jn = 2, nptr                                    ! by sub-basins
106               z3d(nx,:,:) =  ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) 
107               DO jk = 2, jpkm1 
108                  z3d(nx,:,jk) = z3d(nx,:,jk-1) + z3d(nx,:,jk)    ! effective j-Stream-Function (MSF)
109               END DO
110               cl1 = TRIM('zomsf'//clsubb(jn) )
111               CALL iom_put( cl1, z3d * rc_sv )
112            END DO
113         ENDIF
114         !
115      ELSE
116         !
117         IF( iom_use("zotemglo") ) THEN    ! i-mean i-k-surface
118            DO jk = 1, jpkm1
119               DO jj = 1, jpj
120                  DO ji = 1, jpi
121                     zsfc = e1t(ji,jj) * fse3t(ji,jj,jk)
122                     zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
123                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc
124                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc
125                  ENDDO
126               ENDDO
127            ENDDO
128            DO jn = 1, nptr
129               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) )
130               cl1 = TRIM('zosrf'//clsubb(jn) )
131               CALL iom_put( cl1, zmask )
132               !
133               z3d(nx,:,:) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) &
134                  &            / MAX( zmask(1,:,:), 10.e-15 )
135               cl1 = TRIM('zotem'//clsubb(jn) )
136               CALL iom_put( cl1, z3d )
137               !
138               z3d(nx,:,:) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) &
139                  &            / MAX( zmask(1,:,:), 10.e-15 )
140               cl1 = TRIM('zosal'//clsubb(jn) )
141               CALL iom_put( cl1, z3d )
142            END DO
143         ENDIF
144         !
145         !                                ! Advective and diffusive heat and salt transport
146         IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN   
147            z2d(nx,:) = htr_adv(:) * rc_pwatt        !  (conversion in PW)
148            cl1 = 'sophtadv'                 
149            CALL iom_put( TRIM(cl1), z2d )
150            z2d(nx,:) = str_adv(:) * rc_ggram        ! (conversion in Gg)
151            cl1 = 'sopstadv'
152            CALL iom_put( TRIM(cl1), z2d )
153         ENDIF
154         !
155         IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN   
156            z2d(nx,:) = htr_ldf(:) * rc_pwatt        !  (conversion in PW)
157            cl1 = 'sophtldf'
158            CALL iom_put( TRIM(cl1), z2d )
159            z2d(nx,:) = str_ldf(:) * rc_ggram        !  (conversion in Gg)
160            cl1 = 'sopstldf'
161            CALL iom_put( TRIM(cl1), z2d )
162         ENDIF
163         !
164      ENDIF
165      !
166      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr')
167      !
168   END SUBROUTINE dia_ptr
169
170
171   SUBROUTINE dia_ptr_init
172      !!----------------------------------------------------------------------
173      !!                  ***  ROUTINE dia_ptr_init  ***
174      !!                   
175      !! ** Purpose :   Initialization, namelist read
176      !!----------------------------------------------------------------------
177      INTEGER :: jn,  inum   ! local integers
178      !!----------------------------------------------------------------------
179
180      l_diaptr = iom_use("zomsfglo") .OR. iom_use("zotemglo") .OR. iom_use("sophtadv") .OR. iom_use("sophtldf")
181      CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  )
182      IF( inum > 0 ) l_subbas  = .TRUE.
183     
184      IF(lwp) THEN                     ! Control print
185         WRITE(numout,*)
186         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
187         WRITE(numout,*) '~~~~~~~~~~~~'
188         WRITE(numout,*) '   Namelist namptr : set ptr parameters'
189         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr
190         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      l_subbas  = ', l_subbas
191      ENDIF
192
193      IF( l_diaptr ) THEN 
194         !
195         CALL dom_ngb( 180., 90., nx, ny, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots)
196         !
197         IF( l_subbas ) THEN
198            nptr = 5            ! Global, Atlantic, Pacific, Indian, Indo-Pacific
199            ALLOCATE( clsubb(nptr) )
200            clsubb(1) = 'glo' ;  clsubb(2) = 'atl'  ;  clsubb(3) = 'pac'  ;  clsubb(4) = 'ind'  ;  clsubb(5) = 'ipc'
201         ELSE               
202            nptr = 1       ! Global only
203            ALLOCATE( clsubb(nptr) )
204            clsubb(1) = 'glo' 
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( l_subbas ) THEN                ! load sub-basin mask
215            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
216            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
217            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
218            CALL iom_close( inum )
219            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
220            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
221            ELSE WHERE                     ;   btm30(:,:) = ssmask(:,:)
222            END WHERE
223         ENDIF
224   
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         ! Initialise arrays to zero because diatpr is called before they are first calculated
232         ! Note that this means diagnostics will not be exactly correct when model run is restarted.
233         htr_adv(:) = 0._wp  ;  str_adv(:) =  0._wp 
234         htr_ldf(:) = 0._wp  ;  str_ldf(:) =  0._wp 
235         !
236      ENDIF 
237      !
238   END SUBROUTINE dia_ptr_init
239
240
241   FUNCTION dia_ptr_alloc()
242      !!----------------------------------------------------------------------
243      !!                    ***  ROUTINE dia_ptr_alloc  ***
244      !!----------------------------------------------------------------------
245      INTEGER               ::   dia_ptr_alloc   ! return value
246      INTEGER, DIMENSION(3) ::   ierr
247      !!----------------------------------------------------------------------
248      ierr(:) = 0
249      !
250      ALLOCATE( btmsk(jpi,jpj,nptr) ,           &
251         &      htr_adv(jpj) , str_adv(jpj) ,   &
252         &      htr_ldf(jpj) , str_ldf(jpj) , STAT=ierr(1)  )
253         !
254      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2))
255      !
256      ALLOCATE( btm30(jpi,jpj), STAT=ierr(3)  )
257
258         !
259      dia_ptr_alloc = MAXVAL( ierr )
260      IF(lk_mpp)   CALL mpp_sum( dia_ptr_alloc )
261      !
262   END FUNCTION dia_ptr_alloc
263
264
265   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval )
266      !!----------------------------------------------------------------------
267      !!                    ***  ROUTINE ptr_sj_3d  ***
268      !!
269      !! ** Purpose :   i-k sum computation of a j-flux array
270      !!
271      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
272      !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
273      !!
274      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
275      !!----------------------------------------------------------------------
276      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pva   ! mask flux array at V-point
277      REAL(wp), INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
278      !
279      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
280      INTEGER                  ::   ijpj         ! ???
281      REAL(wp), POINTER, DIMENSION(:) :: p_fval  ! function value
282      !!--------------------------------------------------------------------
283      !
284      p_fval => p_fval1d
285
286      ijpj = jpj
287      p_fval(:) = 0._wp
288      IF( PRESENT( pmsk ) ) THEN
289         DO jk = 1, jpkm1
290            DO jj = 2, jpjm1
291               DO ji = fs_2, fs_jpim1   ! Vector opt.
292                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) * pmsk(ji,jj)
293               END DO
294            END DO
295         END DO
296      ELSE
297         DO jk = 1, jpkm1
298            DO jj = 2, jpjm1
299               DO ji = fs_2, fs_jpim1   ! Vector opt.
300                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) 
301               END DO
302            END DO
303         END DO
304      ENDIF
305#if defined key_mpp_mpi
306      IF(lk_mpp)   CALL mpp_sum( p_fval, ijpj, ncomm_znl)
307#endif
308      !
309   END FUNCTION ptr_sj_3d
310
311
312   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval )
313      !!----------------------------------------------------------------------
314      !!                    ***  ROUTINE ptr_sj_2d  ***
315      !!
316      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array
317      !!
318      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
319      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
320      !!
321      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
322      !!----------------------------------------------------------------------
323      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pva   ! mask flux array at V-point
324      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
325      !
326      INTEGER                  ::   ji,jj       ! dummy loop arguments
327      INTEGER                  ::   ijpj        ! ???
328      REAL(wp), POINTER, DIMENSION(:) :: p_fval ! function value
329      !!--------------------------------------------------------------------
330      !
331      p_fval => p_fval1d
332
333      ijpj = jpj
334      p_fval(:) = 0._wp
335      IF( PRESENT( pmsk ) ) THEN
336         DO jj = 2, jpjm1
337            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
338               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) * pmsk(ji,jj)
339            END DO
340         END DO
341      ELSE
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)
345            END DO
346         END DO
347      ENDIF
348#if defined key_mpp_mpi
349      CALL mpp_sum( p_fval, ijpj, ncomm_znl )
350#endif
351      !
352   END FUNCTION ptr_sj_2d
353
354
355   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
356      !!----------------------------------------------------------------------
357      !!                    ***  ROUTINE ptr_sjk  ***
358      !!
359      !! ** Purpose :   i-sum computation of an array
360      !!
361      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
362      !!
363      !! ** Action  : - p_fval: i-mean poleward flux of pva
364      !!----------------------------------------------------------------------
365      !!
366      IMPLICIT none
367      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pta    ! mask flux array at V-point
368      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask
369      !!
370      INTEGER                           :: ji, jj, jk ! dummy loop arguments
371      REAL(wp), POINTER, DIMENSION(:,:) :: p_fval     ! return function value
372#if defined key_mpp_mpi
373      INTEGER, DIMENSION(1) ::   ish
374      INTEGER, DIMENSION(2) ::   ish2
375      INTEGER               ::   ijpjjpk
376      REAL(wp), DIMENSION(jpj*jpk) ::   zwork    ! mask flux array at V-point
377#endif
378      !!--------------------------------------------------------------------
379      !
380      p_fval => p_fval2d
381
382      p_fval(:,:) = 0._wp
383      !
384      IF( PRESENT( pmsk ) ) THEN
385         DO jk = 1, jpkm1
386            DO jj = 2, jpjm1
387!!gm here, use of tmask_i  ==> no need of loop over nldi, nlei....
388               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
389                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj)
390               END DO
391            END DO
392         END DO
393      ELSE
394         DO jk = 1, jpkm1
395            DO jj = 2, jpjm1
396               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
397                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * tmask_i(ji,jj)
398               END DO
399            END DO
400         END DO
401      END IF
402      !
403#if defined key_mpp_mpi
404      ijpjjpk = jpj*jpk
405      ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk
406      zwork(1:ijpjjpk) = RESHAPE( p_fval, ish )
407      CALL mpp_sum( zwork, ijpjjpk, ncomm_znl )
408      p_fval(:,:) = RESHAPE( zwork, ish2 )
409#endif
410      !
411   END FUNCTION ptr_sjk
412
413
414   !!======================================================================
415END MODULE diaptr
Note: See TracBrowser for help on using the repository browser.