source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90 @ 6140

Last change on this file since 6140 was 6140, checked in by timgraham, 5 years ago

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

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