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/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90 @ 9294

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

Merge changes from dev_r5020_CNRS_DIAPTR into trunk, see ticket #1435

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