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.
ptr.F90 in trunk/NEMO/OPA_SRC/DIA – NEMO

source: trunk/NEMO/OPA_SRC/DIA/ptr.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.1 KB
Line 
1MODULE ptr
2   !!======================================================================
3   !!                       ***  MODULE  ptr  ***
4   !! Ocean physics:  brief description of the purpose of the module
5   !!                 (please no more than 2 lines)
6   !!=====================================================================
7#if   defined key_diaptr
8   !!----------------------------------------------------------------------
9   !!   'key_diaptr'  :                    DIAgnostics: Poleward TRansports
10   !!----------------------------------------------------------------------
11   !!   exa_mpl      : liste of module subroutine (caution, never use the
12   !!   exa_mpl_init : name of the module for a routine)
13   !!   exa_mpl_stp  : Please try to use 3 letter block for routine names
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE oce           ! ocean dynamics and active tracers
17   USE dom_oce       ! ocean space and time domain
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! *  Routine accessibility
23   PUBLIC dia_ptr    ! call by stp routine
24   PUBLIC prt_vj     ! call by tra_ldf & tra_adv routines
25
26   !! * Share Module variables
27   LOGICAL, PUBLIC, PARAMETER ::   lk_diaptr = .TRUE.    ! poleward transport flag
28   INTEGER, PUBLIC ::      !!! ** ptr namelist (namptr) **
29      nf_ptr = 15           ! frequency of ptr computation
30   REAL(wp), PUBLIC, DIMENSION(jpj) ::   &   ! poleward transport
31      pht_adv, pst_adv,  &  ! heat and salt: advection
32      pht_ove, pst_ove,  &  ! heat and salt: overturning
33      pht_ldf, pst_ldf,  &  ! heat and salt: lateral diffusion
34      pht_eiv, pst_eiv      ! heat and salt: bolus advection
35
36   !! Module variables
37   REAL(wp), DIMENSION(jpj,jpk) ::   & 
38      tn_jk  , sn_jk  ,  &  ! "zonal" mean temperature and salinity
39      v_msf           ,  &  ! "meridional" Stream-Function
40#if defined key_diaeiv
41      v_msf_eiv       ,  &  ! bolus "meridional" Stream-Function
42#endif
43      surf_jk_r             ! inverse of the ocean "zonal" section surface
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47   !!----------------------------------------------------------------------
48   !!   OPA 9.0 , LODYC-IPSL  (2003)
49   !!----------------------------------------------------------------------
50
51CONTAINS
52
53   FUNCTION ptr_vj( pva )   RESULT ( p_fval )
54      !!----------------------------------------------------------------------
55      !!                    ***  ROUTINE ptr_vj  ***
56      !!
57      !! ** Purpose :   "zonal" and vertical sum computation of a "meridional"
58      !!      flux array
59      !!
60      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
61      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
62      !!
63      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
64      !!
65      !! History :
66      !!   9.0  !  03-09  (G. Madec)  Original code
67      !!----------------------------------------------------------------------
68      !! * arguments
69      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   &
70         pva                         ! mask flux array at V-point
71
72      !! * local declarations
73      INTEGER  ::   ji, jj, jk        ! dummy loop arguments
74      REAL(wp),DIMENSION(jpj) ::   &
75         p_fval                       ! function value
76      !!--------------------------------------------------------------------
77
78      p_fval(  1  ) = 0.e0
79      p_fval(jpjm1) = 0.e0
80      DO jk = 1, jpkm1
81         DO jj = 2, jpjm1
82            DO ji = fs_2, fs_jpim1   ! Vector opt.
83               p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) 
84            END DO
85         END DO
86      END DO
87
88#if defined key_mpp
89      CALL mpp_sum( p_fval, jpj )     !!bug  I presume
90#endif
91
92   END FUNCTION ptr_vj
93
94
95   FUNCTION ptr_vjk( pva )   RESULT ( p_fval )
96      !!----------------------------------------------------------------------
97      !!                    ***  ROUTINE ptr_vjk  ***
98      !!
99      !! ** Purpose :   "zonal" sum computation of a "meridional" flux array
100      !!
101      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
102      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
103      !!
104      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
105      !!
106      !! History :
107      !!   9.0  !  03-09  (G. Madec)  Original code
108      !!----------------------------------------------------------------------
109      !! * arguments
110      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   &
111         pva                         ! mask flux array at V-point
112
113      !! * local declarations
114      INTEGER  ::   ji, jj, jk        ! dummy loop arguments
115      REAL(wp),DIMENSION(jpj,jpk) ::   &
116         p_fval                       ! return function value
117      !!--------------------------------------------------------------------
118 
119      p_fval(  1  , : ) = 0.e0
120      p_fval(jpjm1, : ) = 0.e0
121      p_fval(  :  ,jpk) = 0.e0
122      DO jk = 1, jpkm1
123         DO jj = 2, jpjm1
124            DO ji = fs_2, fs_jpim1   ! Vector opt.
125               p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj)
126            END DO
127         END DO
128      END DO
129
130#if defined key_mpp
131      CALL mpp_sum( p_fval, jpj*jpk )    !!bug  I presume
132#endif
133
134   END FUNCTION ptr_vjk
135
136   FUNCTION ptr_vtjk( pva )   RESULT ( p_fval )
137      !!----------------------------------------------------------------------
138      !!                    ***  ROUTINE ptr_vtjk  ***
139      !!
140      !! ** Purpose :   "zonal" mean computation of a tracer field
141      !!
142      !! ** Method  : - i-sum of mj(pva) using the interior 2D vmask (vmask_i)
143      !!      multiplied by the inverse of the surface of the "zonal" ocean
144      !!      section
145      !!
146      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
147      !!
148      !! History :
149      !!   9.0  !  03-09  (G. Madec)  Original code
150      !!----------------------------------------------------------------------
151      !! * arguments
152      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   &
153         pva                         ! mask flux array at V-point
154 
155      !! * local declarations
156      INTEGER  ::   ji, jj, jk        ! dummy loop arguments
157      REAL(wp),DIMENSION(jpj,jpk) ::   &
158         p_fval                       ! return function value
159      !!--------------------------------------------------------------------
160      p_fval(  1  , : ) = 0.e0
161      p_fval(jpjm1, : ) = 0.e0
162      p_fval(  :  ,jpk) = 0.e0
163      DO jk = 1, jpkm1
164         DO jj = 2, jpjm1
165            DO ji = fs_2, fs_jpim1   ! Vector opt.
166               p_fval(jj,jk) = p_fval(jj,jk) + ( pva(ji,jj,jk) + pva(ji,jj,jk) )                &
167                  &                          * e1v(ji,jj) * fse3v(ji,jj,jk) * vmask(ji,jj,jk)   &
168                  &                          * tmask_i(ji,jj+1) * tmask_i(ji,jj)
169            END DO
170         END DO
171      END DO
172      p_fval(:,:) = p_val(:,:) * 0.5
173
174#if defined key_mpp
175      CALL mpp_sum( p_fval, jpj*jpk )         !!bug  I presume
176#endif
177
178   END FUNCTION ptr_vtjk
179
180
181   SUBROUTINE dia_ptr( kt )
182      !!----------------------------------------------------------------------
183      !!                  ***  ROUTINE dia_ptr  ***
184      !!----------------------------------------------------------------------
185      !! * Argument
186      INTEGER, INTENT(in) ::   kt   ! ocean time step index
187
188      !! * Local variables
189      REAL(wp) ::   &
190         zsverdrup = 1.e-6          ! conversion from m3/s to Sverdrup
191      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
192         zvfl                       ! mask flux array at V-point
193      !!----------------------------------------------------------------------
194
195      ! "zonal" mean temperature and salinity at V-points
196      tn_jk(:,:) = prt_vtjk( tn(:,:,:) ) * surf_jk_r(:,:)
197      sn_jk(:,:) = prt_vtjk( sn(:,:,:) ) * surf_jk_r(:,:)
198
199      ! "zonal" mean mass flux at V-points
200      v_msf(:,:) = prt_vjk( vn(:,:,:) ) 
201#if defined key_diaeiv
202      ! "zonal" mean bolus mass flux at V-points
203      v_msf_eiv(:,:) = prt_vjk( v_eiv(:,:,:) ) 
204      ! Bolus "Meridional" Stream-Function
205      DO jk = jpkm1, 1
206         v_msf_eiv(:,jk) = v_msf_eiv(:,jk-1) + v_msf_eiv(:,jk)
207      END DO
208      v_msf_eiv(:,:) = v_msf_eiv(:,:) * zsverdrup
209#endif
210
211      ! poleward transport: overturning component
212      pht_ove(:) = SUM( v_msf(:,:) * tn_jk(:,:), 2 )   ! SUM over jk
213
214      ! "Meridional" Stream-Function
215      DO jk = jpkm1, 1
216         v_msf(:,jk) = v_msf(:,jk-1) + v_msf(:,jk)
217      END DO
218      v_msf(:,:) = v_msf(:,:) * zsverdrup
219
220      ! output
221      CALL dia_ptr_wri( kt )
222
223   END SUBROUTINE dia_ptr
224
225
226   SUBROUTINE dia_ptr_init
227      !!----------------------------------------------------------------------
228      !!                  ***  ROUTINE dia_ptr_init  ***
229      !!                   
230      !! ** Purpose :   initialization of ....
231      !!
232      !! ** Method  :   blah blah blah ...
233      !!
234      !! ** input   :   Namlist namptr
235      !!
236      !! ** Action  :   ... 
237      !!
238      !! history :
239      !!   9.0  !  03-08  (Autor Names)  Original code
240      !!----------------------------------------------------------------------
241      !! * local declarations
242      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
243         z_1             ! temporary workspace
244
245      NAMELIST/namptr/ nf_ptr
246      !!----------------------------------------------------------------------
247
248      ! Read Namelist namptr : poleward transport parameters
249      REWIND ( numptr )
250      READ   ( numnam, namptr )
251
252
253      ! Control print
254      IF(lwp) THEN
255         WRITE(numout,*)
256         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization'
257         WRITE(numout,*) '~~~~~~~~~~~~'
258         WRITE(numout,*) '          Namelist namptr : set ptr parameters'
259         WRITE(numout,*) '             frequency of computation       nf_ptr  = ', nf_ptr
260      ENDIF
261
262      ! inverse of the ocean "zonal" v-point section
263      z_1(:,:,:) = 1.e0
264      surf_jk_r(:,:) = prt_vtjk( z1(:,:,:) )
265      WHERE( surf_jk_r(:,:) /= 0.e0 )   surf_jk_r(:,:) = 1.e0 / surf_jk_r(:,:)
266
267   END SUBROUTINE dia_ptr_init
268
269#if defined key_fdir
270   !!---------------------------------------------------------------------
271   !!   'key_fdir'                                      direct access file
272   !!---------------------------------------------------------------------
273
274   SUBROUTINE dia_ptr_wri( kt )
275      !!---------------------------------------------------------------------
276      !!                ***  ROUTINE diaptr  ***
277      !!
278      !! ** Purpose :   output of poleward fluxes
279      !!
280      !! ** Method  :   NetCDF file
281      !!
282      !! History :
283      !!   9.0  !  03-09  (G. Madec)  Original code 
284      !!----------------------------------------------------------------------
285      !! * Arguments
286      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
287      !!----------------------------------------------------------------------
288
289      IF( kt == nit000 ) THEN
290
291         ! Reference latitude
292         ! ------------------
293         !                                           ! =======================
294         IF( cp_cfg == "orca" ) THEN                 !   ORCA configurations
295            !                                        ! =======================
296 
297            iline = 100 / jp_cfg             ! i-line that passes near the North Pole
298            zphi(:) = 0.e0   
299            DO ji = mi0(iline), mi1(iline)
300               zphi(:) = gphiv(ji,:)         ! if iline is in the local domain
301            END DO
302#  if defined key_mpp
303            CALL mpp_sum( zphi, jpj )        ! provide the correct zphi to all local domains
304#  endif
305            !                                        ! =======================
306         ELSE                                        !   OTHER configurations
307            !                                        ! =======================
308            zphi(:) = gphiv(1,:)             ! assume lat/lon coordinate, select the first i-line
309            !
310         ENDIF
311 
312         ! open the output file         
313         CALL ctlopn( numptr, 'opaptr.output', 'UNKNOWN', 'UNFORMATTED', 'SEQUENTIAL', 1, numout, lw
314p, 1 )
315
316         ! header of output
317         WRITE( numptr ) cexper, no, zdt, nf_ptr, jpj, jpk, zphi
318
319      ENDIF
320
321      IF( MOD( kt, nf_ptr ) == 0 ) THEN
322            IF(lwp) WRITE( numptr ) kt, tn_jk, sn_jk, v_msf,              &
323# if defined key_diaeiv
324               &                    pht_eiv, pst_eiv, v_msf_eiv,          &
325# endif
326               &                    pht_adv, pht_ldf, pst_adv, pst_ldf
327      ENDIF
328
329      IF( kt === nitend )   CLOSE( numptr ) 
330
331   END SUBROUTINE dia_ptr_wri
332
333#else
334   !!---------------------------------------------------------------------
335   !!   Default option :                                       NetCDF file
336   !!---------------------------------------------------------------------
337
338   SUBROUTINE dia_ptr_wri( kt )
339      !!---------------------------------------------------------------------
340      !!                ***  ROUTINE diaptr  ***
341      !!
342      !! ** Purpose :   output of poleward fluxes
343      !!
344      !! ** Method  :   NetCDF file
345      !!
346      !! History :
347      !!   9.0  !  03-09  (G. Madec)  Original code
348      !!----------------------------------------------------------------------
349      USE ioipsl          ! NetCDF IPSL library
350
351      !! * Arguments
352      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
353
354      !! * Save variables   
355      INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw   &
356          , ndex(1)
357
358      !! * Local variables
359      CHARACTER (len=15) ::   clexp
360      CHARACTER (len=40) ::   &
361         clhstnam, clop, clmax      ! temporary names
362      REAL(wp) ::   &
363         zsto, zout, zdt,           ! temporary scalars
364         zjulian
365      !!----------------------------------------------------------------------
366     
367      ! Define frequency of output and means
368      zdt = rdt
369      IF( nacc == 1 ) zdt = rdtmin
370#if defined key_diainstant
371         zsto = nf_ptr * zdt
372         clop = "inst(x)"               ! no use of the mask value (require less cpu time)
373         !!! clop="inst(only(x))"       ! put 1.e+20 on land (very expensive!!)
374#else
375         zsto = zdt
376         clop = "ave(x)"                ! no use of the mask value (require less cpu time)
377         !!! clop="ave(only(x))"        ! put 1.e+20 on land (very expensive!!)
378#endif
379      zout = nf_ptr * zdt
380      zmax = ( nitend - nit000 + 1 ) * zdt
381     
382         
383      ! define time axis
384      it = kt - nit000 + 1
385
386      ! Initialization
387      ! --------------
388      IF( kt == nit000 ) THEN
389     
390      zdt = rdt
391      IF( nacc == 1 ) zdt = rdtmin
392
393         ! Reference latitude
394         ! ------------------
395         !                                           ! =======================
396         IF( cp_cfg == "orca" ) THEN                 !   ORCA configurations
397            !                                        ! =======================
398
399            iline = 100 / jp_cfg             ! i-line that passes near the North Pole
400            zphi(:) = 0.e0                   
401            DO ji = mi0(iline), mi1(iline) 
402               zphi(:) = gphiv(ji,:)         ! if iline is in the local domain
403            END DO
404#  if defined key_mpp
405            CALL mpp_sum( zphi, jpj )        ! provide the correct zphi to all local domains
406#  endif
407            !                                        ! =======================
408         ELSE                                        !   OTHER configurations
409            !                                        ! =======================
410            zphi(:) = gphiv(1,:)             ! assume lat/lon coordinate, select the first i-line
411            !
412         ENDIF
413
414         ! OPEN netcdf file
415         ! ----------------
416         ! Define frequency of output and means
417         zsto = nf_ptr * zdt
418         clop = "ave(x)"
419         zout = nf_ptr * zdt
420         zfoo = 0.e0
421
422         ! Compute julian date from starting date of the run
423
424         CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian )
425
426         CALL dia_nam( clhstnam, nf_ptr, 'diaptr' )
427         IF(lwp)WRITE( numout,*)" Name of diaptr NETCDF file ",clhstnam
428
429         ! Horizontal grid : zphi()
430         CALL histbeg(clhstnam, 1, zfoo, jpj, zphi,   &
431            1, 1, 1, jpj, 0, zjulian, zdt, nhoridz, numptr )
432         ! Vertical grids : gdept, gdepw
433         CALL histvert( numptr, "deptht", "Vertical T levels",   &
434            "m", jpk, gdept, ndepidzt )
435         CALL histvert( numptr, "depthw", "Vertical W levels",   &
436            "m", jpk, gdepw, ndepidzw )
437         
438         !  Zonal mean T and S
439         
440         CALL histdef( numptr, "zotemglo", "Zonal Mean Temperature","C" ,   &
441            1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
442         CALL histdef( numptr, "zosalglo", "Zonal Mean Salinity","PSU"  ,   &
443            1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout )
444
445         !  Meridional Stream-Function (eulerian and bolus)
446         
447         CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: global","Sv" ,   &
448            1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
449
450         !  Heat transport
451
452         CALL histdef( numptr, "sophtadv", "Advective Heat Transport"      ,   &
453            "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
454         CALL histdef( numptr, "sophtldf", "Diffusive Heat Transport"      ,   &
455            "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
456         CALL histdef( numptr, "sophtove", "Overturning Heat Transport"    ,   &
457            "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
458
459         !  Salt transport
460
461         CALL histdef( numptr, "sopstadv", "Advective Salt Transport"      ,   &
462            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
463         CALL histdef( numptr, "sopstldf", "Diffusive Salt Transport"      ,   &
464            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
465         CALL histdef( numptr, "sopstove", "Overturning Salt Transport"    ,   &
466            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
467
468#if defined key_diaeiv
469         ! Eddy induced velocity
470         CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: global",   &
471            "Sv"      , 1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout )
472         CALL histdef( numptr, "sophteiv", "Bolus Advective Heat Transport",   &
473            "PW"      , 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
474         CALL histdef( numptr, "sopsteiv", "Bolus Advective Salt Transport",   &
475            "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout )
476#endif
477
478         CALL histend( numptr )
479
480      ENDIF
481
482      IF( MOD( kt, nf_ptr ) == 0 ) THEN
483
484         ! define time axis
485         it= kt - nit000 + 1
486         CALL histwrite( numptr, "zotemglo", it, tn_jk    , jpj*jpk, ndex )
487         CALL histwrite( numptr, "zosalglo", it, sn_jk    , jpj*jpk, ndex )
488         CALL histwrite( numptr, "zomsfglo", it, v_msf    , jpj*jpk, ndex )
489         CALL histwrite( numptr, "sophtadv", it, pht_adv  , jpj    , ndex )
490         CALL histwrite( numptr, "sophtldf", it, pht_ldf  , jpj    , ndex )
491         CALL histwrite( numptr, "sophtove", it, pht_ove  , jpj    , ndex )
492         CALL histwrite( numptr, "sopstadv", it, pst_adv  , jpj    , ndex )
493         CALL histwrite( numptr, "sopstldf", it, pst_ldf  , jpj    , ndex )
494         CALL histwrite( numptr, "sopstove", it, pst_ove  , jpj    , ndex )
495#if defined key_diaeiv
496         CALL histwrite( numptr, "zomsfeiv", it, v_msf_eiv, jpj*jpk, ndex )
497         CALL histwrite( numptr, "sophteiv", it, pht_eiv  , jpj    , ndex )
498         CALL histwrite( numptr, "sopsteiv", it, pst_eiv  , jpj    , ndex )
499#endif
500 
501      ENDIF
502
503      ! Close the file
504      IF( kt === nitend )   CALL histclo( numptr )           ! Netcdf write
505
506   END SUBROUTINE dia_ptr_wri
507
508#endif
509
510#else
511   !!----------------------------------------------------------------------
512   !!   Default option :                                       Empty Module
513   !!----------------------------------------------------------------------
514   LOGICAL, PUBLIC, PARAMETER ::   lk_diaptr = .FALSE.    ! poleward transport flag
515CONTAINS
516   SUBROUTINE dia_ptr( kt )        ! Empty routine
517      WRITE(*,*) kt
518   END SUBROUTINE dia_ptr
519#endif
520
521   !!======================================================================
522END MODULE ptr
Note: See TracBrowser for help on using the repository browser.