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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

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