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 @ 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: 17.5 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_fsc && ! defined key_coupled
7   !!----------------------------------------------------------------------
8   !!   "key_dynspg_fsc" 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 , LODYC-IPSL  (2003)
46   !!----------------------------------------------------------------------
47
48CONTAINS
49
50   SUBROUTINE dia_fwb( kt )
51      !!---------------------------------------------------------------------
52      !!                  ***  ROUTINE dia_fwb  ***
53      !!     
54      !! ** Purpose :
55      !!
56      !! ** Method :
57      !!
58      !! History :
59      !!   8.2  !  01-02  (E. Durand)  Original code
60      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
61      !!----------------------------------------------------------------------
62      !! * Arguments
63      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
64
65      !! * Local declarations
66      INTEGER :: ji, jj, jk, jt   ! dummy loop indices
67      REAL(wp) ::   zarea, zvol, zwei
68      REAL(wp) ::  ztemi(4), ztemo(4), zsali(4), zsalo(4), zflxi(4), zflxo(4)
69      REAL(wp) ::  zt, zs, zu 
70      REAL(wp) ::  zsm0, zempnew
71      !!----------------------------------------------------------------------
72
73      ! Mean global salinity
74      zsm0 = 34.72654
75
76      ! To compute emp mean value mean emp
77
78      IF( kt == nit000 ) THEN
79
80         a_emp    = 0.e0
81         a_precip = 0.e0
82         a_rnf    = 0.e0
83         a_sshb   = 0.e0 ! valeur de ssh au debut de la simulation
84         a_salb   = 0.e0 ! valeur de sal au debut de la simulation
85         a_aminus = 0.e0
86         a_aplus  = 0.e0
87         ! sshb used because diafwb called after tranxt (i.e. after the swap)
88         a_sshb = SUM( e1t(:,:) * e2t(:,:) * sshb(:,:) * tmask_i(:,:) )
89
90         DO jk = 1, jpkm1
91            DO jj = 2, jpjm1
92               DO ji = fs_2, fs_jpim1   ! vector opt.
93                  zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
94                  a_salb = a_salb + ( sb(ji,jj,jk) - zsm0 ) * zwei
95               END DO
96            END DO
97         END DO
98      ENDIF
99     
100      a_emp    = SUM( e1t(:,:) * e2t(:,:) * emp   (:,:) * tmask_i(:,:) )
101#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
102      a_precip = SUM( e1t(:,:) * e2t(:,:) * watm  (:,:) * tmask_i(:,:) )
103#endif
104      a_rnf    = SUM( e1t(:,:) * e2t(:,:) * runoff(:,:) * tmask_i(:,:) )
105
106      IF( aminus /= 0.0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
107      IF( aplus  /= 0.0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
108
109#if defined key_mpp
110      ! Mpp: sum over all the global domain
111      CALL  mpp_sum( a_sshn )                            !!!!!! bugggggg   a_sshn note befined before!!!!!
112#endif
113
114      IF( kt == nitend ) THEN
115         a_sshn = 0.e0
116         a_saln = 0.e0
117         zarea = 0.e0
118         zvol  = 0.e0
119         zempnew = 0.e0
120         ! Mean sea level at nitend
121         a_sshn = SUM( e1t(:,:) * e2t(:,:) * sshn(:,:) * tmask_i(:,:) )
122         zarea  = SUM( e1t(:,:) * e2t(:,:) *             tmask_i(:,:) )
123         
124         DO jk = 1, jpkm1   
125            DO jj = 2, jpjm1
126               DO ji = fs_2, fs_jpim1   ! vector opt.
127                  zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
128                  a_saln = a_saln + ( sn(ji,jj,jk) - zsm0 ) * zwei
129                  zvol  = zvol  + zwei
130               END DO
131            END DO
132         END DO
133         
134         a_aminus = a_aminus/(nitend-nit000+1)
135         a_aplus  = a_aplus/(nitend-nit000+1)
136
137         ! Conversion in m3
138         a_emp    = a_emp * rdttra(1) * 1.e-3 
139         a_precip = a_precip * rdttra(1) * 1.e-3 / rday
140         a_rnf    = a_rnf * rdttra(1) * 1.e-3
141         
142         ! Alpha1=Alpha0-Rest/(Precip+runoff)
143         !  C A U T I O N : precipitations are negative !!     
144         
145         zempnew = a_sshn / ( ( nitend - nit000 + 1 ) * rdt ) * 1.e3 / zarea
146
147      ENDIF
148
149
150      ! Calcul des termes de transport
151      ! ------------------------------
152     
153      ! 1 --> Gibraltar
154      ! 2 --> Cadiz
155      ! 3 --> Red Sea
156      ! 4 --> Baltic Sea
157
158      IF( kt == nit000 ) THEN
159         a_flxi(:) = 0.e0
160         a_flxo(:) = 0.e0
161         a_temi(:) = 0.e0
162         a_temo(:) = 0.e0
163         a_sali(:) = 0.e0
164         a_salo(:) = 0.e0
165      ENDIF
166
167      zflxi(:) = 0.e0
168      zflxo(:) = 0.e0
169      ztemi(:) = 0.e0
170      ztemo(:) = 0.e0
171      zsali(:) = 0.e0
172      zsalo(:) = 0.e0
173
174      ! Mean flow at Gibraltar
175
176      IF( cp_cfg == "orca" ) THEN 
177               
178         SELECT CASE ( jp_cfg )
179         !                                           ! =======================
180         CASE ( 4 )                                  !  ORCA_R4 configuration
181            !                                        ! =======================
182            ji =  mi1(70)
183            jj =  mj1(52)
184            !                                        ! =======================
185         CASE ( 2 )                                  !  ORCA_R2 configuration
186            !                                        ! =======================
187            ji = mi1(139)
188            jj = mj1(102)
189            !                                        ! =======================
190         CASE DEFAULT                                !    ORCA R05 or R025
191            !                                        ! =======================
192            IF(lwp) WRITE(numout,cform_err)
193            IF(lwp) WRITE(numout,*)' dia_fwb Not yet implemented in ORCA_R05 or R025'
194            nstop = nstop + 1
195            !
196         END SELECT
197         !
198      ENDIF
199      DO jk = 1, 18 
200         zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
201         zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
202         zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
203
204         IF( un(ji,jj,jk) > 0.0 ) THEN
205            zflxi(1) = zflxi(1) +    zu
206            ztemi(1) = ztemi(1) + zt*zu
207            zsali(1) = zsali(1) + zs*zu
208         ELSE
209            zflxo(1) = zflxo(1) +    zu
210            ztemo(1) = ztemo(1) + zt*zu
211            zsalo(1) = zsalo(1) + zs*zu
212         ENDIF
213      END DO
214     
215      ! Mean flow at Cadiz
216      IF( cp_cfg == "orca" ) THEN
217               
218         SELECT CASE ( jp_cfg )
219         !                                           ! =======================
220         CASE ( 4 )                                  !  ORCA_R4 configuration
221            !                                        ! =======================
222            ji =  mi1(69  )
223            jj =  mj1(52  )
224            !                                        ! =======================
225         CASE ( 2 )                                  !  ORCA_R2 configuration
226            !                                        ! =======================
227            ji = mi1(137)
228            jj = mj1(102)
229            !                                        ! =======================
230         CASE DEFAULT                                !    ORCA R05 or R025
231            !                                        ! =======================
232            IF(lwp) WRITE(numout,cform_err)
233            IF(lwp) WRITE(numout,*)' dia_fwb Not yet implemented in ORCA_R05 or R025'
234            nstop = nstop + 1
235            !
236         END SELECT
237         !
238      ENDIF
239      DO jk = 1, 23 
240         zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
241         zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
242         zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
243         
244         IF( un(ji,jj,jk) > 0.0 ) THEN
245            zflxi(2) = zflxi(2) +    zu
246            ztemi(2) = ztemi(2) + zt*zu
247            zsali(2) = zsali(2) + zs*zu
248         ELSE
249            zflxo(2) = zflxo(2) +    zu
250            ztemo(2) = ztemo(2) + zt*zu
251            zsalo(2) = zsalo(2) + zs*zu
252         ENDIF
253      END DO
254
255      ! Mean flow at Red Sea entrance
256      IF( cp_cfg == "orca" ) THEN
257               
258         SELECT CASE ( jp_cfg )
259         !                                           ! =======================
260         CASE ( 4 )                                  !  ORCA_R4 configuration
261            !                                        ! =======================
262            ji =  mi1(83  )
263            jj =  mj1(45  )
264            !                                        ! =======================
265         CASE ( 2 )                                  !  ORCA_R2 configuration
266            !                                        ! =======================
267            ji = mi1(161)
268            jj = mj1( 88)
269            !                                        ! =======================
270         CASE DEFAULT                                !    ORCA R05 or R025
271            !                                        ! =======================
272            IF(lwp) WRITE(numout,cform_err)
273            IF(lwp) WRITE(numout,*)' dia_fwb Not yet implemented in ORCA_R05 or R025'
274            nstop = nstop + 1
275            !
276         END SELECT
277         !
278      ENDIF
279      DO jk = 1, 15 
280         zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
281         zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
282         zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
283         
284         IF( un(ji,jj,jk) > 0.0 ) THEN
285            zflxi(3) = zflxi(3) +    zu
286            ztemi(3) = ztemi(3) + zt*zu
287            zsali(3) = zsali(3) + zs*zu
288         ELSE
289            zflxo(3) = zflxo(3) +    zu
290            ztemo(3) = ztemo(3) + zt*zu
291            zsalo(3) = zsalo(3) + zs*zu
292         ENDIF
293      END DO
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            ji =   1     ! Not in the domain
303            jj =   1 
304            !                                        ! =======================
305         CASE ( 2 )                                  !  ORCA_R2 configuration
306            !                                        ! =======================
307            ji = mi1(146)
308            jj = mj1(116)
309            !                                        ! =======================
310         CASE DEFAULT                                !    ORCA R05 or R025
311            !                                        ! =======================
312            IF(lwp) WRITE(numout,cform_err)
313            IF(lwp) WRITE(numout,*)' dia_fwb Not yet implemented in ORCA_R05 or R025'
314            nstop = nstop + 1
315            !
316         END SELECT
317         !
318      ENDIF
319      DO jk = 1, 20
320         zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
321         zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
322         zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
323         
324         IF( un(ji,jj,jk) > 0.0 ) THEN
325            zflxi(4) = zflxi(4) +    zu
326            ztemi(4) = ztemi(4) + zt*zu
327            zsali(4) = zsali(4) + zs*zu
328         ELSE
329            zflxo(4) = zflxo(4) +    zu
330            ztemo(4) = ztemo(4) + zt*zu
331            zsalo(4) = zsalo(4) + zs*zu
332         ENDIF
333      END DO
334
335      ! Sum at each time-step
336      DO jt = 1, 4 
337         IF( zflxi(jt) /= 0.0 .AND. zflxo(jt) /= 0.0 ) THEN
338            a_flxi(jt) = a_flxi(jt) + zflxi(jt)
339            a_temi(jt) = a_temi(jt) + ztemi(jt)/zflxi(jt)
340            a_sali(jt) = a_sali(jt) + zsali(jt)/zflxi(jt)
341            a_flxo(jt) = a_flxo(jt) + zflxo(jt)
342            a_temo(jt) = a_temo(jt) + ztemo(jt)/zflxo(jt)
343            a_salo(jt) = a_salo(jt) + zsalo(jt)/zflxo(jt)
344         ENDIF
345      END DO
346
347      IF( kt == nitend ) THEN
348         DO jt = 1, 4 
349            a_flxi(jt) = a_flxi(jt)/((nitend-nit000+1)*1.e6)
350            a_temi(jt) = a_temi(jt)/( nitend-nit000+1)
351            a_sali(jt) = a_sali(jt)/( nitend-nit000+1)
352            a_flxo(jt) = a_flxo(jt)/((nitend-nit000+1)*1.e6)
353            a_temo(jt) = a_temo(jt)/( nitend-nit000+1)
354            a_salo(jt) = a_salo(jt)/( nitend-nit000+1)
355         END DO
356      ENDIF
357
358
359      ! Ecriture des diagnostiques
360      ! --------------------------
361
362      IF ( kt == nitend ) THEN
363
364         OPEN(111,FILE='STRAIT.dat')
365         WRITE(111,*)
366         WRITE(111,*)    'Net freshwater budget '
367         WRITE(111,9010) '  emp    = ',a_emp,   ' m3 =', a_emp   /((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
368         WRITE(111,9010) '  precip = ',a_precip,' m3 =', a_precip/((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
369         WRITE(111,9010) '  a_rnf   = ',a_rnf,   ' m3 =', a_rnf   /((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
370         WRITE(111,*)
371         WRITE(111,9010) '  zarea =',zarea
372         WRITE(111,9010) '  zvol  =',zvol
373         WRITE(111,*)
374         WRITE(111,*)    'Mean sea level : '
375         WRITE(111,9010) '  at nit000 = ',a_sshb        ,' m3 '
376         WRITE(111,9010) '  at nitend = ',a_sshn        ,' m3 '
377         WRITE(111,9010) '  diff      = ',(a_sshn-a_sshb),' m3 =', (a_sshn-a_sshb)/((nitend-nit000+1)*rdt) * 1.e-6,' Sv'
378         WRITE(111,9020) '  mean sea level elevation    =', a_sshn/zarea,' m'
379         WRITE(111,*)
380         WRITE(111,*)    'Anomaly of salinity content : '
381         WRITE(111,9010) '  at nit000 = ',a_salb        ,' psu.m3 '
382         WRITE(111,9010) '  at nitend = ',a_saln        ,' psu.m3 '
383         WRITE(111,9010) '  diff      = ',(a_saln-a_salb),' psu.m3'
384         WRITE(111,*)
385         WRITE(111,*)    'Mean salinity : '
386         WRITE(111,9020) '  at nit000 =',a_salb/zvol+zsm0   ,' psu '
387         WRITE(111,9020) '  at nitend =',a_saln/zvol+zsm0   ,' psu '
388         WRITE(111,9020) '  diff      =',(a_saln-a_salb)/zvol,' psu'
389         WRITE(111,9020) '  S-SLevitus=',a_saln/zvol,' psu'
390         WRITE(111,*)
391         WRITE(111,*)    'Coeff : '
392         WRITE(111,9030) '  Alpha+   =  ', a_aplus
393         WRITE(111,9030) '  Alpha-   =  ', a_aminus
394         WRITE(111,*)
395         WRITE(111,*)
396         WRITE(111,*)    'Gibraltar : '
397         WRITE(111,9030) '  Flux entrant (Sv) :', a_flxi(1)
398         WRITE(111,9030) '  Flux sortant (Sv) :', a_flxo(1)
399         WRITE(111,9030) '  T entrant (deg)   :', a_temi(1)
400         WRITE(111,9030) '  T sortant (deg)   :', a_temo(1)
401         WRITE(111,9030) '  S entrant (psu)   :', a_sali(1)
402         WRITE(111,9030) '  S sortant (psu)   :', a_salo(1)
403         WRITE(111,*)
404         WRITE(111,*)    'Cadiz : '
405         WRITE(111,9030) '  Flux entrant (Sv) :', a_flxi(2)
406         WRITE(111,9030) '  Flux sortant (Sv) :', a_flxo(2)
407         WRITE(111,9030) '  T entrant (deg)   :', a_temi(2)
408         WRITE(111,9030) '  T sortant (deg)   :', a_temo(2)
409         WRITE(111,9030) '  S entrant (psu)   :', a_sali(2)
410         WRITE(111,9030) '  S sortant (psu)   :', a_salo(2)
411         WRITE(111,*)
412         WRITE(111,*)    'Bab el Mandeb : '
413         WRITE(111,9030) '  Flux entrant (Sv) :', a_flxi(3)
414         WRITE(111,9030) '  Flux sortant (Sv) :', a_flxo(3)
415         WRITE(111,9030) '  T entrant (deg)   :', a_temi(3)
416         WRITE(111,9030) '  T sortant (deg)   :', a_temo(3)
417         WRITE(111,9030) '  S entrant (psu)   :', a_sali(3)
418         WRITE(111,9030) '  S sortant (psu)   :', a_salo(3)
419         WRITE(111,*)
420         WRITE(111,*)    'Baltic : '
421         WRITE(111,9030) '  Flux entrant (Sv) :', a_flxi(4)
422         WRITE(111,9030) '  Flux sortant (Sv) :', a_flxo(4)
423         WRITE(111,9030) '  T entrant (deg)   :', a_temi(4)
424         WRITE(111,9030) '  T sortant (deg)   :', a_temo(4)
425         WRITE(111,9030) '  S entrant (psu)   :', a_sali(4)
426         WRITE(111,9030) '  S sortant (psu)   :', a_salo(4)
427         CLOSE(111)
428      ENDIF
429
430 9005 FORMAT(1X,A,ES24.16)
431 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
432 9020 FORMAT(1X,A,F10.5,A)
433 9030 FORMAT(1X,A,F8.2,A)
434
435   END SUBROUTINE dia_fwb
436
437#else
438   !!----------------------------------------------------------------------
439   !!   Default option :                                       Empty Module
440   !!----------------------------------------------------------------------
441   LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .FALSE.    ! fresh water budget flag
442CONTAINS
443   SUBROUTINE dia_fwb( kt )        ! Empty routine
444      WRITE(*,*) kt                      ! no warning in compilation phase
445   END SUBROUTINE dia_fwb
446#endif
447
448   !!======================================================================
449END MODULE diafwb
Note: See TracBrowser for help on using the repository browser.