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/NEMO/OPA_SRC/DIA – NEMO

source: trunk/NEMO/OPA_SRC/DIA/diafwb.F90 @ 405

Last change on this file since 405 was 389, checked in by opalod, 18 years ago

RB:nemo_v1_update_038: first integration of Agrif :

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