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 branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90 @ 3973

Last change on this file since 3973 was 3973, checked in by clevy, 11 years ago

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