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.
diafwb_gm.F90 on Ticket #127 – Attachment – NEMO

Ticket #127: diafwb_gm.F90

File diafwb_gm.F90, 19.1 KB (added by gm, 15 years ago)

bug correction for mpp from Gurvan

Line 
1MODULE diafwb
2   !!======================================================================
3   !!                       ***  MODULE  diafwb  ***
4   !! Ocean diagnostics: freshwater budget
5   !!======================================================================
6   !! History :  OPA  !  2001-02  (E. Durand)  Original code
7   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
8   !!            3.2  !  2009-07  (A.R. Porter, G. Madec) correct mpp bugs
9   !!----------------------------------------------------------------------
10#if ( defined key_orca_r2 || defined  key_orca_r4 ) && ! defined key_dynspg_rl && ! defined key_coupled
11   !!----------------------------------------------------------------------
12   !!   NOT "key_dynspg_rl" and "key_orca_r2 or 4"
13   !!----------------------------------------------------------------------
14   !!----------------------------------------------------------------------
15   !!   dia_fwb     : freshwater budget for global ocean configurations
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE phycst          ! physical constants
20   USE sbc_oce         ! ???
21   USE zdf_oce         ! ocean vertical physics
22   USE in_out_manager  ! I/O manager
23   USE lib_mpp         ! distributed memory computing library
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC dia_fwb    ! routine called by step.F90
29
30   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .TRUE.    !: fresh water budget flag
31
32   REAL(wp)               ::   a_emp ,          &
33      &                        a_sshb, a_sshn, a_salb, a_saln
34   REAL(wp), DIMENSION(4) ::   a_flxi, a_flxo, a_temi, a_temo, a_sali, a_salo
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !!   OPA 9.0 , LOCEAN-IPSL (2006)
41   !! $Id$
42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE dia_fwb( kt )
48      !!---------------------------------------------------------------------
49      !!                  ***  ROUTINE dia_fwb  ***
50      !!     
51      !! ** Purpose :
52      !!----------------------------------------------------------------------
53      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
54      !!
55      INTEGER :: inum             ! temporary logical unit
56      INTEGER :: ji, jj, jk, jt   ! dummy loop indices
57      INTEGER :: ii0, ii1, ij0, ij1
58      REAL(wp) ::   zarea, zvol, zwei
59      REAL(wp) ::  ztemi(4), ztemo(4), zsali(4), zsalo(4), zflxi(4), zflxo(4)
60      REAL(wp) ::  zt, zs, zu 
61      REAL(wp) ::  zsm0
62      !!----------------------------------------------------------------------
63
64      ! Mean global salinity
65      zsm0 = 34.72654
66
67      ! To compute emp mean value mean emp
68
69      IF( kt == nit000 ) THEN
70
71         a_emp    = 0.e0
72         a_sshb   = 0.e0 ! valeur de ssh au debut de la simulation
73         a_salb   = 0.e0 ! valeur de sal au debut de la simulation
74         ! sshb used because diafwb called after tranxt (i.e. after the swap)
75         a_sshb = SUM( e1t(:,:) * e2t(:,:) * sshb(:,:) * tmask_i(:,:) )
76         IF( lk_mpp )   CALL mpp_sum( a_sshb )      ! sum over the global domain
77
78         DO jk = 1, jpkm1
79            DO jj = 2, jpjm1
80               DO ji = fs_2, fs_jpim1   ! vector opt.
81                  zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
82                  a_salb = a_salb + ( sb(ji,jj,jk) - zsm0 ) * zwei
83               END DO
84            END DO
85         END DO
86         IF( .NOT. lk_vvl )   a_salb = a_salb + SUM( e1t(:,:) * e2t(:,:) * sshb(:,:) * tmask_i(:,:) * ( sb(:,:,1) - zsm0 ) )
87         IF( lk_mpp )   CALL mpp_sum( a_salb )      ! sum over the global domain
88      ENDIF
89     
90      a_emp    = SUM( e1t(:,:) * e2t(:,:) * emp   (:,:) * tmask_i(:,:) )
91      IF( lk_mpp )   CALL mpp_sum( a_emp    )       ! sum over the global domain
92
93      IF( kt == nitend ) THEN
94         a_sshn = 0.e0
95         a_saln = 0.e0
96         zarea = 0.e0
97         zvol  = 0.e0
98         ! Mean sea level at nitend
99         a_sshn = SUM( e1t(:,:) * e2t(:,:) * sshn(:,:) * tmask_i(:,:) )
100         IF( lk_mpp )   CALL mpp_sum( a_sshn )      ! sum over the global domain
101         zarea  = SUM( e1t(:,:) * e2t(:,:) *             tmask_i(:,:) )
102         IF( lk_mpp )   CALL mpp_sum( zarea  )      ! sum over the global domain
103         
104         DO jk = 1, jpkm1   
105            DO jj = 2, jpjm1
106               DO ji = fs_2, fs_jpim1   ! vector opt.
107                  zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
108                  a_saln = a_saln + ( sn(ji,jj,jk) - zsm0 ) * zwei
109                  zvol  = zvol  + zwei
110               END DO
111            END DO
112         END DO
113         IF( .NOT. lk_vvl ) THEN                    ! correct with salt content in the ssh 'layer'
114            a_saln = a_saln + SUM(  e1t(ji,jj) * e2t(ji,jj) * tmask_i(:,:) * ssh(:,:) * ( sn(:,:,1) - zsm0 )  )
115            zvol   = zvol   + SUM(  e1t(ji,jj) * e2t(ji,jj) * tmask_i(:,:) * ssh(:,:)                         )
116         ENDIF
117         IF( lk_mpp )   CALL mpp_sum( a_saln )      ! sum over the global domain
118         IF( lk_mpp )   CALL mpp_sum( zvol )        ! sum over the global domain
119         
120         ! Conversion in m3
121         a_emp    = a_emp * rdttra(1) * 1.e-3 
122         !
123      ENDIF
124
125
126      ! Calcul des termes de transport
127      ! ------------------------------
128     
129      ! 1 --> Gibraltar
130      ! 2 --> Cadiz
131      ! 3 --> Red Sea
132      ! 4 --> Baltic Sea
133
134      IF( kt == nit000 ) THEN
135         a_flxi(:) = 0.e0
136         a_flxo(:) = 0.e0
137         a_temi(:) = 0.e0
138         a_temo(:) = 0.e0
139         a_sali(:) = 0.e0
140         a_salo(:) = 0.e0
141      ENDIF
142
143      zflxi(:) = 0.e0
144      zflxo(:) = 0.e0
145      ztemi(:) = 0.e0
146      ztemo(:) = 0.e0
147      zsali(:) = 0.e0
148      zsalo(:) = 0.e0
149
150      ! Mean flow at Gibraltar
151
152      IF( cp_cfg == "orca" ) THEN 
153               
154         SELECT CASE ( jp_cfg )
155         !                                           ! =======================
156         CASE ( 4 )                                  !  ORCA_R4 configuration
157            !                                        ! =======================
158            ii0 = 70   ;   ii1 = 70
159            ij0 = 52   ;   ij1 = 52
160            !                                        ! =======================
161         CASE ( 2 )                                  !  ORCA_R2 configuration
162            !                                        ! =======================
163            ii0 = 140   ;   ii1 = 140
164            ij0 = 102   ;   ij1 = 102
165            !                                        ! =======================
166         CASE DEFAULT                                !    ORCA R05 or R025
167            !                                        ! =======================
168            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
169            !
170         END SELECT
171         !
172         DO jj = mj0(ij0), mj1(ij1)
173            DO ji = mi0(ii0), MAX( mi1(ii1), jpim1 )      ! max used to avoid mpp out-of-bound
174               DO jk = 1, jpk 
175                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
176                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
177                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)   ! only if interior point
178
179                  IF( un(ji,jj,jk) > 0.e0 ) THEN
180                     zflxi(1) = zflxi(1) +    zu
181                     ztemi(1) = ztemi(1) + zt*zu
182                     zsali(1) = zsali(1) + zs*zu
183                  ELSE
184                     zflxo(1) = zflxo(1) +    zu
185                     ztemo(1) = ztemo(1) + zt*zu
186                     zsalo(1) = zsalo(1) + zs*zu
187                  ENDIF
188               END DO
189            END DO
190         END DO
191      ENDIF
192     
193      ! Mean flow at Cadiz
194      IF( cp_cfg == "orca" ) THEN
195               
196         SELECT CASE ( jp_cfg )
197         !                                           ! =======================
198         CASE ( 4 )                                  !  ORCA_R4 configuration
199            !                                        ! =======================
200            ii0 = 69   ;   ii1 = 69
201            ij0 = 52   ;   ij1 = 52
202            !                                        ! =======================
203         CASE ( 2 )                                  !  ORCA_R2 configuration
204            !                                        ! =======================
205            ii0 = 137   ;   ii1 = 137
206            ij0 = 101   ;   ij1 = 102
207            !                                        ! =======================
208         CASE DEFAULT                                !    ORCA R05 or R025
209            !                                        ! =======================
210            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
211            !
212         END SELECT
213         !
214         DO jj = mj0(ij0), mj1(ij1)
215            DO ji = mi0(ii0), MAX( mi1(ii1), jpim1 )      ! max used to avoid mpp out-of-bound
216               DO jk = 1, jpk 
217                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
218                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
219                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)   ! only if interior point
220                 
221                  IF( un(ji,jj,jk) > 0.e0 ) THEN
222                     zflxi(2) = zflxi(2) +    zu
223                     ztemi(2) = ztemi(2) + zt*zu
224                     zsali(2) = zsali(2) + zs*zu
225                  ELSE
226                     zflxo(2) = zflxo(2) +    zu
227                     ztemo(2) = ztemo(2) + zt*zu
228                     zsalo(2) = zsalo(2) + zs*zu
229                  ENDIF
230               END DO
231            END DO
232         END DO
233      ENDIF
234
235      ! Mean flow at Red Sea entrance
236      IF( cp_cfg == "orca" ) THEN
237               
238         SELECT CASE ( jp_cfg )
239         !                                           ! =======================
240         CASE ( 4 )                                  !  ORCA_R4 configuration
241            !                                        ! =======================
242            ii0 = 83   ;   ii1 = 83
243            ij0 = 45   ;   ij1 = 45
244            !                                        ! =======================
245         CASE ( 2 )                                  !  ORCA_R2 configuration
246            !                                        ! =======================
247            ii0 = 160   ;   ii1 = 160
248            ij0 = 88    ;   ij1 = 88 
249            !                                        ! =======================
250         CASE DEFAULT                                !    ORCA R05 or R025
251            !                                        ! =======================
252            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
253            !
254         END SELECT
255         !
256         DO jj = mj0(ij0), mj1(ij1)
257            DO ji = mi0(ii0), MAX( mi1(ii1), jpim1 )      ! max used to avoid mpp out-of-bound
258               DO jk = 1, jpk 
259                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
260                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
261                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)   ! only if interior point
262                 
263                  IF( un(ji,jj,jk) > 0.e0 ) THEN
264                     zflxi(3) = zflxi(3) +    zu
265                     ztemi(3) = ztemi(3) + zt*zu
266                     zsali(3) = zsali(3) + zs*zu
267                  ELSE
268                     zflxo(3) = zflxo(3) +    zu
269                     ztemo(3) = ztemo(3) + zt*zu
270                     zsalo(3) = zsalo(3) + zs*zu
271                  ENDIF
272               END DO
273            END DO
274         END DO
275      ENDIF
276
277      ! Mean flow at Baltic Sea entrance
278      IF( cp_cfg == "orca" ) THEN
279               
280         SELECT CASE ( jp_cfg )
281         !                                           ! =======================
282         CASE ( 4 )                                  !  ORCA_R4 configuration
283            !                                        ! =======================
284            ii0 = 1     ;   ii1 = 1 
285            ij0 = 1     ;   ij1 = 1 
286            !                                        ! =======================
287         CASE ( 2 )                                  !  ORCA_R2 configuration
288            !                                        ! =======================
289            ii0 = 146   ;   ii1 = 146 
290            ij0 = 116   ;   ij1 = 116
291            !                                        ! =======================
292         CASE DEFAULT                                !    ORCA R05 or R025
293            !                                        ! =======================
294            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
295            !
296         END SELECT
297         !
298         DO jj = mj0(ij0), mj1(ij1)
299            DO ji = mi0(ii0), MAX( mi1(ii1), jpim1 )      ! max used to avoid mpp out-of-bound
300               DO jk = 1, jpk
301                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
302                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
303                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)   ! only if interior point
304                 
305                  IF( un(ji,jj,jk) > 0.e0 ) THEN
306                     zflxi(4) = zflxi(4) +    zu
307                     ztemi(4) = ztemi(4) + zt*zu
308                     zsali(4) = zsali(4) + zs*zu
309                  ELSE
310                     zflxo(4) = zflxo(4) +    zu
311                     ztemo(4) = ztemo(4) + zt*zu
312                     zsalo(4) = zsalo(4) + zs*zu
313                  ENDIF
314               END DO
315            END DO
316         END DO
317      ENDIF
318
319      ! Sum at each time-step
320      DO jt = 1, 4 
321         !
322         IF( zflxi(jt) /= 0.e0 ) THEN
323            a_flxi(jt) = a_flxi(jt) + zflxi(jt)
324            a_temi(jt) = a_temi(jt) + ztemi(jt)/zflxi(jt)
325            a_sali(jt) = a_sali(jt) + zsali(jt)/zflxi(jt)
326         ENDIF
327         !
328         IF( zflxo(jt) /= 0.e0 ) THEN
329            a_flxo(jt) = a_flxo(jt) + zflxo(jt)
330            a_temo(jt) = a_temo(jt) + ztemo(jt)/zflxo(jt)
331            a_salo(jt) = a_salo(jt) + zsalo(jt)/zflxo(jt)
332         ENDIF
333         !
334      END DO
335
336      IF( kt == nitend ) THEN
337         DO jt = 1, 4 
338            a_flxi(jt) = a_flxi(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
339            a_temi(jt) = a_temi(jt) /   FLOAT( nitend - nit000 + 1 )
340            a_sali(jt) = a_sali(jt) /   FLOAT( nitend - nit000 + 1 )
341            a_flxo(jt) = a_flxo(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
342            a_temo(jt) = a_temo(jt) /   FLOAT( nitend - nit000 + 1 )
343            a_salo(jt) = a_salo(jt) /   FLOAT( nitend - nit000 + 1 )
344         END DO
345         IF( lk_mpp ) THEN
346            CALL mpp_sum( a_flxi, 4 )      ! sum over the global domain
347            CALL mpp_sum( a_temi, 4 )      ! sum over the global domain
348            CALL mpp_sum( a_sali, 4 )      ! sum over the global domain
349
350            CALL mpp_sum( a_flxo, 4 )      ! sum over the global domain
351            CALL mpp_sum( a_temo, 4 )      ! sum over the global domain
352            CALL mpp_sum( a_salo, 4 )      ! sum over the global domain
353         ENDIF
354      ENDIF
355
356
357      ! Ecriture des diagnostiques
358      ! --------------------------
359
360      IF ( kt == nitend .AND. cp_cfg == "orca" ) THEN
361
362         CALL ctlopn( inum, 'STRAIT.dat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
363            &         1, numout, lwp, 1 )
364         WRITE(inum,*)
365         WRITE(inum,*)    'Net freshwater budget '
366         WRITE(inum,9010) '  emp    = ',a_emp,   ' m3 =', a_emp   /(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
367         WRITE(inum,*)
368         WRITE(inum,9010) '  zarea =',zarea
369         WRITE(inum,9010) '  zvol  =',zvol
370         WRITE(inum,*)
371         WRITE(inum,*)    'Mean sea level : '
372         WRITE(inum,9010) '  at nit000 = ',a_sshb        ,' m3 '
373         WRITE(inum,9010) '  at nitend = ',a_sshn        ,' m3 '
374         WRITE(inum,9010) '  diff      = ',(a_sshn-a_sshb),' m3 =', (a_sshn-a_sshb)/(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
375         WRITE(inum,9020) '  mean sea level elevation    =', a_sshn/zarea,' m'
376         WRITE(inum,*)
377         WRITE(inum,*)    'Anomaly of salinity content : '
378         WRITE(inum,9010) '  at nit000 = ',a_salb        ,' psu.m3 '
379         WRITE(inum,9010) '  at nitend = ',a_saln        ,' psu.m3 '
380         WRITE(inum,9010) '  diff      = ',(a_saln-a_salb),' psu.m3'
381         WRITE(inum,*)
382         WRITE(inum,*)    'Mean salinity : '
383         WRITE(inum,9020) '  at nit000 =',a_salb/zvol+zsm0   ,' psu '
384         WRITE(inum,9020) '  at nitend =',a_saln/zvol+zsm0   ,' psu '
385         WRITE(inum,9020) '  diff      =',(a_saln-a_salb)/zvol,' psu'
386         WRITE(inum,9020) '  S-SLevitus=',a_saln/zvol,' psu'
387         WRITE(inum,*)
388         WRITE(inum,*)    'Gibraltar : '
389         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(1)
390         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(1)
391         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(1)
392         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(1)
393         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(1)
394         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(1)
395         WRITE(inum,*)
396         WRITE(inum,*)    'Cadiz : '
397         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(2)
398         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(2)
399         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(2)
400         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(2)
401         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(2)
402         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(2)
403         WRITE(inum,*)
404         WRITE(inum,*)    'Bab el Mandeb : '
405         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(3)
406         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(3)
407         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(3)
408         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(3)
409         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(3)
410         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(3)
411         WRITE(inum,*)
412         WRITE(inum,*)    'Baltic : '
413         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(4)
414         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(4)
415         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(4)
416         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(4)
417         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(4)
418         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(4)
419         CLOSE(inum)
420      ENDIF
421
422 9005 FORMAT(1X,A,ES24.16)
423 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
424 9020 FORMAT(1X,A,F10.5,A)
425 9030 FORMAT(1X,A,F9.4,A)
426
427   END SUBROUTINE dia_fwb
428
429#else
430   !!----------------------------------------------------------------------
431   !!   Default option :                                       Dummy Module
432   !!----------------------------------------------------------------------
433   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .FALSE.    !: fresh water budget flag
434CONTAINS
435   SUBROUTINE dia_fwb( kt )        ! Empty routine
436      WRITE(*,*) 'dia_fwb: : You should not have seen this print! error?', kt
437   END SUBROUTINE dia_fwb
438#endif
439
440   !!======================================================================
441END MODULE diafwb