source: trunk/SOURCES/source_3.20/ice_sal_diff_CW.f @ 4

Last change on this file since 4 was 4, checked in by vancop, 8 years ago

initial import /Users/ioulianikolskaia/Boulot/CODES/LIM1D/ARCHIVE/TMP/LIM1D_v3.20/

File size: 15.6 KB
Line 
1      SUBROUTINE ice_sal_diff_CW(nlay_i,kideb,kiut)
2
3      !!------------------------------------------------------------------
4      !!                ***         ROUTINE ice_sal_diff      ***
5      !!
6      !! ** Purpose :
7      !!        This routine computes new salinities in the ice
8      !!
9      !! ** Method  : Vertical salinity profile computation
10      !!              Resolves brine transport equation
11      !!           
12      !! ** Steps
13      !!
14      !! ** Arguments
15      !!
16      !! ** Inputs / Outputs
17      !!
18      !! ** External
19      !!
20      !! ** References : Vancop. et al., 2008
21      !!
22      !! ** History  :
23      !!    (06-2003) Martin Vancop. LIM1D
24      !!    (06-2008) Martin Vancop. BIO-LIM
25      !!    (09-2008) Martin Vancop. Explicit gravity drainage
26      !!
27      !!------------------------------------------------------------------
28
29      USE lib_fortran
30
31      INCLUDE 'type.com'
32      INCLUDE 'para.com'
33      INCLUDE 'const.com'
34      INCLUDE 'ice.com'
35      INCLUDE 'thermo.com'
36
37      REAL(8), DIMENSION(nlay_i) :: 
38     &   z_ms_i             ,    !: mass of salt times thickness
39     &   z_sbr_i                 !: brine salinity
40
41      REAL(8), DIMENSION(nlay_i) ::  !: dummy factors for tracer equation
42     &   za                 ,    !: winter
43     &   zb                 ,   
44     &   ze                 ,    !: summer
45     &   zf                 ,   
46     &   zind               ,    !: independent term in the tridiag system
47     &   zindw              ,    !: independent term in the tridiag system
48     &   zinds              ,    !: independent term in the tridiag system
49     &   zindtbis           ,    !:
50     &   zdiagbis           ,    !:
51     &   zflux                   !: flux of tracer under layer i
52
53      REAL(8), DIMENSION(nlay_i,3) :: !: dummy factors for tracer equation
54     &   ztrid              ,    !: tridiagonal matrix
55     &   ztridw             ,    !: tridiagonal matrix, winter
56     &   ztrids                  !: tridiagonal matrix, summer
57
58      REAL(8) :: 
59     &   zdummy1            ,    !: dummy factors
60     &   zdummy2            ,    !:
61     &   zdummy3            ,    !:
62     &   zswitch_open       ,    !: switch for brine network open or not
63     &   zswitchw           ,    !: switch for winter drainage
64     &   zswitchs           ,    !: switch for summer drainage
65     &   zeps      = 1.0e-20     !: numerical limit
66
67      ! Rayleigh number computation
68      REAL(8) ::
69     &   ze_i_min           ,    !: minimum brine volume
70     &   zc                 ,    !: temporary scalar for sea ice specific heat
71     &   zk                 ,    !: temporary scalar for sea ice thermal conductivity
72     &   zalphara                !: multiplicator for diffusivity
73
74      REAL(8), DIMENSION(nlay_i) :: 
75     &   zsigma             ,    !: brine salinity at layer interfaces
76     &   zperm              ,    !: permeability
77     &   zpermin            ,    !: minimum permeability
78     &   zrhodiff           ,    !: density difference
79     &   zlevel             ,    !: height of the water column
80     &   zthdiff            ,    !: thermal diffusivity
81     &   zgrad_t                 !: temperature gradient in the ice
82
83      INTEGER :: 
84     &   layer2             ,    !: layer loop index
85     &   indtr                   !: index of tridiagonal system
86
87      CHARACTER(len=4)      :: 
88     &   bc = 'conc'             !: Boundary condition 'conc' or 'flux'
89
90      REAL(8) :: 
91     &   z_ms_i_ini         ,    !: initial mass of salt
92     &   z_ms_i_fin         ,    !: final mass of salt
93     &   z_fs_b             ,    !: basal flux of salt
94     &   z_fs_su            ,    !: surface flux of salt
95     &   z_dms_i                 !: mass variation       
96
97      LOGICAL ::
98     &   ln_write           ,
99     &   ln_con             ,
100     &   ln_sal             ,
101     &   ln_be              ,
102     &   ln_gd              ,
103     &   ln_fl
104
105      ln_write = .TRUE.  ! write outputs
106      ln_con   = .TRUE.  ! conservation check
107      ln_sal   = .TRUE.  ! compute salinity variations or not
108      ln_be    = .TRUE.  ! compute brine expulsion or not
109      ln_gd    = .TRUE.  ! compute gravity drainage or not
110      ln_fl    = .TRUE.  ! compute flushing or not
111
112      IF ( ln_write ) THEN
113         WRITE(numout,*)
114         WRITE(numout,*) ' ** ice_sal_diff_CW : '
115         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ '
116         WRITE(numout,*) ' Cox and weeks based gravity drainage '
117         WRITE(numout,*) ' ln_sal = ', ln_sal
118      ENDIF
119
120      IF ( ln_sal ) THEN
121!
122!------------------------------------------------------------------------------|
123! 1) Initialization
124!------------------------------------------------------------------------------|
125!
126      IF ( ln_write ) THEN
127         WRITE(numout,*) ' - Initialization ... '
128      ENDIF
129
130      DO 10 ji = kideb, kiut
131
132      ! gravity drainage parameter (Cox and Weeks 88)
133      zeta            = 20.
134
135      ! brine diffusivity
136      diff_br(:) = 0.0
137
138      !---------------------------
139      ! Brine volume and salinity
140      !---------------------------
141      DO layer = 1, nlay_i
142         e_i_b(layer) = - tmut * s_i_b(ji,layer) / ( t_i_b(ji,layer) 
143     &                - tpw )
144         z_sbr_i(layer) = s_i_b(ji,layer) / e_i_b(layer)
145      END DO
146
147      !----------
148      ! Switches
149      !----------
150      ! summer switch = 1 if Tsu ge tpw and min brine volume superior than e_thr_flu
151      zswitchs = MAX( 0.0, SIGN ( 1.0d0, t_su_b(ji) - tpw ) ) ! 0 si hiver 1 si ete
152      zbvmin   = 1.0
153      DO layer = 1, nlay_i
154         zbvmin = MIN( e_i_b(layer) , zbvmin ) ! minimum brine volume
155      END DO
156      IF ( zbvmin .LT. e_thr_flu ) zswitchs = 0.0
157
158      ! winter switch
159      zswitchw = 1.0 - zswitchs
160
161      !------------------
162      ! Percolating flux
163      !------------------
164      ! Percolating flow ( rho dh * beta * switch / rhow )
165      qsummer = ( - rhog * MIN ( dh_i_surf(ji) , 0.0 ) 
166     &            - rhon * MIN ( dh_s_tot(ji)  , 0.0 ) ) 
167      qsummer = qsummer * flu_beta * zswitchs / 1000.0
168
169      !--------------------
170      ! Conservation check
171      !--------------------
172      IF ( ln_con ) THEN
173         CALL ice_sal_column( kideb , kiut , z_ms_i_ini , 
174     &                    s_i_b(1,1:nlay_i),
175     &                    deltaz_i_phy, nlay_i, .FALSE. )
176      ENDIF ! ln_con
177
178      IF ( ln_write ) THEN
179         WRITE(numout,*) ' nlay_i      : ', nlay_i
180         WRITE(numout,*) ' kideb       : ', kideb
181         WRITE(numout,*) ' kiut        : ', kiut 
182         WRITE(numout,*) ' '
183         WRITE(numout,*) ' deltaz_i_phy : ', ( deltaz_i_phy(layer), 
184     &                   layer = 1, nlay_i )
185         WRITE(numout,*) ' z_i_phy      : ', ( z_i_phy(layer), 
186     &                   layer = 1, nlay_i )
187         WRITE(numout,*) ' s_i_b       : ', ( s_i_b    (ji,layer), 
188     &                   layer = 1, nlay_i )
189         WRITE(numout,*) ' t_i_b       : ', ( t_i_b    (ji,layer), 
190     &                   layer = 1, nlay_i )
191         WRITE(numout,*) ' e_i_b       : ', ( e_i_b      (layer), 
192     &                   layer = 1, nlay_i )
193         WRITE(numout,*) ' z_sbr_i     : ', ( z_sbr_i    (layer), 
194     &                   layer = 1, nlay_i )
195         WRITE(numout,*)
196         WRITE(numout,*) ' zswitchs : ', zswitchs
197         WRITE(numout,*) ' zswitchw : ', zswitchw
198         WRITE(numout,*) 
199      ENDIF ! ln_write
200
201 10   CONTINUE
202
203!
204!------------------------------------------------------------------------------|
205! 2) Gravity drainage as from Cox and Weeks (1988)
206!------------------------------------------------------------------------------|
207!
208      IF ( zswitchw .EQ. 1.0 ) THEN
209
210      DO 20 ji = kideb, kiut
211
212         !----------------------
213         ! temperature gradient
214         !----------------------
215         z_h_lay    = ht_i_b(ji) / nlay_i
216         zgrad_t(1) = 2. * ( t_i_b(ji,1) - ( ht_s_b(ji)*t_i_b(ji,1) +
217     &                z_h_lay*t_s_b(ji,1) ) / 
218     &                ( z_h_lay + ht_s_b(ji) ) ) / z_h_lay
219
220         DO layer = 2, nlay_i - 1
221            zgrad_t(layer) = ( t_i_b(ji,layer+1) - t_i_b(ji,layer-1) ) /
222     &                       z_h_lay
223         END DO
224         zgrad_t(nlay_i) = - 2.*( t_i_b(ji,nlay_i) - t_bo_b(ji) ) /
225     &                       z_h_lay
226
227         IF ( ln_write ) WRITE(numout,*) ' zgrad_t : ', 
228     &                   ( zgrad_t(layer), layer = 1, nlay_i )
229
230         igrd = 1 ! switch for gravity drainage ( 1 if yes )
231         zdummysb = - FLOAT(igrd) * rhog / 1000. * ht_i_b(ji) ! salt flux [ kg NaCl.m-2.s-1 ]
232         zdsdt = 0.0
233         DO layer = nlay_i, 1, -1
234            IF ( e_i_b(layer) .LE. e_thr_flu ) igrd = 0
235            IF ( zgrad_t(layer) .LE. 0 )   igrd = 0 ! temperature gradient must
236                                                    ! be directed downwards
237            zds_grd = MIN ( 0.0, delta_cw * ( 1.0 - zeta * e_i_b(layer)
238     &              * zgrad_t(layer) * ddtb * igrd ) )
239            sn_i_b(layer) = s_i_b(ji,layer) + zds_grd
240            zdsdt = zdsdt + FLOAT(igrd) * zds_grd / ddtb / FLOAT(nlay_i)
241         END DO
242         fsb = zdummysb * zdsdt
243
244 20   CONTINUE
245
246      ENDIF ! zswitchw
247
248!
249!------------------------------------------------------------------------------|
250! 3) Compute dummy factors for tracer diffusion equation
251!------------------------------------------------------------------------------|
252
253      IF ( ln_write ) THEN
254         WRITE(numout,*) ' - Compute dummy factors for tracer diffusion'
255         WRITE(numout,*) ' '
256      ENDIF
257
258      DO 30 ji = kideb, kiut 
259         
260      !----------------
261      ! Winter factors
262      !----------------
263      IF ( zswitchw .EQ. 1. ) THEN
264         za(:) = 0.0
265         zb(:) = 0.0
266      ENDIF
267
268      !----------------------
269      ! Summer factors
270      !----------------------
271      ! ze factors
272      DO layer = 1, nlay_i
273         ze(layer) = qsummer * zswitchs / 
274     &               ( e_i_b(layer) * deltaz_i_phy(layer) )
275      END DO ! layer
276
277      ! zf factors
278      ! could remove those, they are totally useless!!! ;-)
279      DO layer = 1, nlay_i - 1
280         zf(layer) = 1./2. * deltaz_i_phy(layer) / 
281     &               ( z_i_phy(layer+1)  - z_i_phy(layer) )
282      END DO ! layer
283
284      IF ( ln_write ) THEN
285         WRITE(numout,*)
286         WRITE(numout,*) ' -Summer factors '
287         WRITE(numout,*) ' ze : ', ( ze(layer), layer = 1, nlay_i ) 
288         WRITE(numout,*) ' zf : ', ( zf(layer), layer = 1, nlay_i ) 
289      ENDIF
290
291 30   CONTINUE
292!
293!-----------------------------------------------------------------------
294! 4) Tridiagonal system terms for tracer diffusion equation, winter
295!-----------------------------------------------------------------------
296!
297      DO 40 ji = kideb, kiut 
298
299         ztridw(:,:) = 0.
300
301 40   CONTINUE
302!
303!-----------------------------------------------------------------------
304! 5) Tridiagonal system terms for tracer diffusion equation, summer
305!-----------------------------------------------------------------------
306!
307      DO 50 ji = kideb, kiut 
308
309      DO layer = 1, nlay_i
310         ztrids(layer,1) = - ze(layer)
311         ztrids(layer,2) = 1.0 + ze(layer)
312         ztrids(layer,3) = 0.0
313         zinds(layer) = z_sbr_i(layer)
314      END DO
315      ztrids(1,1) = 0.0
316
317      IF ( ln_write ) THEN
318         WRITE(numout,*) 
319         WRITE(numout,*) ' -Tridiag terms, summer ... '
320         WRITE(numout,*)
321         DO layer = 1, nlay_i
322            WRITE(numout,*) ' layer : ', layer
323            WRITE(numout,*) ' ztrids   : ', ztrids(layer,1), 
324     &                ztrids(layer,2), ztrids(layer,3)
325            WRITE(numout,*) ' zinds     : ',zinds(layer)
326         END DO
327      ENDIF
328
329 50   CONTINUE
330     
331!
332!-----------------------------------------------------------------------
333! 6) Partitionning tridiag system between summer and winter
334!-----------------------------------------------------------------------
335!
336      DO 60 ji = kideb, kiut 
337
338      DO indtr = 1, 3
339         DO layer = 1, nlay_i
340            ztrid(layer,indtr) = zswitchs * ztrids(layer,indtr)
341         END DO ! layer
342      END DO ! indtr
343
344      DO layer = 1, nlay_i
345         zind(layer) = zswitchs * zinds(layer)
346      END DO ! layer
347
348      IF ( ln_write ) THEN
349         WRITE(numout,*) 
350         WRITE(numout,*) ' -Tridiag terms...'
351         WRITE(numout,*)
352         WRITE(numout,*) ' zswitchw : ', zswitchw
353         WRITE(numout,*) ' zswitchs : ', zswitchs
354         DO layer = 1, nlay_i
355            WRITE(numout,*) ' layer    : ', layer
356            WRITE(numout,*) ' ztrid    : ', ztrid(layer,1), 
357     &                        ztrid(layer,2), ztrid(layer,3)
358            WRITE(numout,*) ' zind     : ', zind(layer)
359         END DO ! layer
360      ENDIF
361
362 60   CONTINUE
363!
364!-----------------------------------------------------------------------
365! 7) Solving the tridiagonal system
366!-----------------------------------------------------------------------
367!
368      DO 70 ji = kideb, kiut 
369
370      ! The tridiagonal system is solved with Gauss elimination
371      ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
372      ! McGraw-Hill 1984.       
373
374      zindtbis(1) =  zind(1)
375      zdiagbis(1) =  ztrid(1,2)
376
377      DO layer = 2, nlay_i
378         zdiagbis(layer)  =  ztrid(layer,2) - ztrid(layer,1) *
379     &                       ztrid(layer-1,3) / zdiagbis(layer-1)
380         zindtbis(layer)  =  zind(layer) - ztrid(layer,1) *
381     &                       zindtbis(layer-1) / zdiagbis(layer-1)
382      END DO
383
384      ! Recover brine salinity
385      z_sbr_i(nlay_i)     =  zindtbis(nlay_i) / zdiagbis(nlay_i) 
386      DO layer = nlay_i - 1 , 1 , -1
387         z_sbr_i(layer)   =  ( zindtbis(layer) - ztrid(layer,3)*
388     &                       z_sbr_i(layer+1)) / zdiagbis(layer)
389      END DO
390      ! Recover ice salinity
391      IF ( zswitchs .EQ. 1.0 ) THEN
392      DO layer = 1, nlay_i
393         sn_i_b(layer)  =  z_sbr_i(layer) * e_i_b(layer)
394!        s_i_b(ji,layer)  =  z_sbr_i(layer) * e_i_b(layer)
395      END DO
396      ENDIF
397
398      IF ( ln_write ) THEN
399         WRITE(numout,*) 
400         WRITE(numout,*) ' -Solving the tridiagonal system ... '
401         WRITE(numout,*)
402         WRITE(numout,*) ' zdiagbis: ', ( zdiagbis(layer) , 
403     &                   layer = 1, nlay_i )
404         WRITE(numout,*) ' zindtbis: ', ( zdiagbis(layer) , 
405     &                   layer = 1, nlay_i )
406         WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i(layer) , 
407     &                   layer = 1, nlay_i )
408         WRITE(numout,*) ' sn_i_b   : ', ( sn_i_b(layer) , 
409     &                   layer = 1, nlay_i )
410      ENDIF
411
412 70   CONTINUE
413!
414!-----------------------------------------------------------------------
415! 8) Mass of salt conserved ?
416!-----------------------------------------------------------------------
417!
418      DO 80 ji = kideb, kiut 
419
420      ! Final mass of salt
421      CALL ice_sal_column( kideb , kiut , z_ms_i_fin , 
422     &                    sn_i_b(1:nlay_i),
423     &                    deltaz_i_phy, nlay_i, .FALSE. )
424
425      ! Bottom flux ( positive upwards )
426      zswitch_open = 0.0
427      IF ( e_i_b(nlay_i) .GE. e_thr_flu ) zswitch_open = 1.0
428      zfb      = zswitchw * ( - e_i_b( nlay_i ) ! had a minus before
429     &                * diff_br(nlay_i) * 2.0
430     &                / deltaz_i_phy(nlay_i) * ( z_sbr_i(nlay_i)
431     &                  - oce_sal ) ) * zswitch_open
432     &                + zswitchs * ( - qsummer * z_sbr_i(nlay_i) )
433     &                / ddtb
434
435      zflux(nlay_i) = zfb
436
437      ! Surface flux of salt
438      zfsu     = zswitchw * 0.0
439
440      ! conservation check
441      CALL ice_sal_conserv(kideb,kiut,'ice_sal_diff : ',zerror,
442     &                           z_ms_i_ini,z_ms_i_fin,
443     &                           zfb   , zfsu  , ddtb)
444
445 80   CONTINUE
446
447      ENDIF ! ln_sal
448!
449!------------------------------------------------------------------------------|
450! End of la sous-routine
451      WRITE(numout,*)
452      END SUBROUTINE 
Note: See TracBrowser for help on using the repository browser.