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.
sedadv.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/TOP_SRC/SED – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/TOP_SRC/SED/sedadv.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 13 years ago

set proper svn properties to all files...

  • Property svn:keywords set to Id
File size: 16.9 KB
Line 
1MODULE sedadv
2#if defined key_sed
3   !!======================================================================
4   !!              ***  MODULE  sedadv  ***
5   !!    Sediment : vertical advection and burial
6   !!=====================================================================
7   !! * Modules used
8   !!----------------------------------------------------------------------
9   !!   sed_adv :
10   !!----------------------------------------------------------------------
11   USE sed     ! sediment global variable
12
13   PUBLIC sed_adv
14
15   !! * Module variable
16   INTEGER, PARAMETER  ::   &
17      nztime = jpksed         ! number of time step between sunrise and sunset
18   !     
19   
20   REAL(wp), DIMENSION(jpksed), SAVE :: &
21      dvolsp, dvolsm, &
22      c2por, ckpor
23
24   REAL(wp), SAVE :: &
25      por1clay , &
26      cpor   
27
28   REAL(wp) :: &
29      eps = 1.e-13
30
31CONTAINS
32
33   SUBROUTINE sed_adv( kt )
34      !!-------------------------------------------------------------------------
35      !!                  ***  ROUTINE sed_adv  ***
36      !!
37      !! ** Purpose : vertical solid sediment advection and burial
38      !!
39      !! ** Method  : At each grid point the 1-dimensional solid sediment column
40      !!              is shifted according the rain added to the top layer and
41      !!              the gaps produced through redissolution so that in the end
42      !!              the original sediment mixed layer geometry is reestablished.
43      !!
44      !!
45      !!   History :
46      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
47      !!        !  04-10 (N. Emprin, M. Gehlen ) F90
48      !!        !  06-04 (C. Ethe)  Re-organization
49      !!-------------------------------------------------------------------------
50      !!* Arguments
51      INTEGER, INTENT(in) ::  &
52         kt                     ! time step
53      ! * local variables
54
55      INTEGER ::  &
56         ji, jk, js
57
58      INTEGER ::  &
59         jn, ntimes, ikwneg
60     
61      REAL(wp), DIMENSION(jpksed,jpsol) :: zsolcpno
62      REAL(wp), DIMENSION(jpksed)       :: zfilled, zfull, zfromup, zempty
63     
64      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zgap, zwb
65      REAL(wp), DIMENSION(jpoce,jpsol)      :: zrainrf
66
67      REAL(wp), DIMENSION(nztime) ::  zraipush
68
69      REAL(wp) :: &
70         zkwnup, zkwnlo, zfrac, &
71         zfromce, zrest
72
73
74      !------------------------------------------------------------------------
75
76
77      IF( kt == nitsed000 ) THEN
78         WRITE(numsed,*) ' '
79         WRITE(numsed,*) ' sed_adv : vertical sediment advection  '
80         WRITE(numsed,*) ' '
81         por1clay = dens * por1(jpksed) * dz(jpksed) / mol_wgt(jsclay)
82         cpor     = por1(jpksed) / por1(2)
83         DO jk = 2, jpksed
84            c2por(jk) = por1(2)      / por1(jk)
85            ckpor(jk) = por1(jpksed) / por1(jk)
86         ENDDO
87         DO jk = jpksedm1, 2, -1
88            dvolsp(jk) = vols(jk+1) / vols(jk)
89         ENDDO
90         DO jk = 3, jpksed
91           dvolsm(jk) = vols(jk-1) / vols(jk)
92        ENDDO
93      ENDIF
94
95
96      ALLOCATE( zgap (jpoce,jpksed) )    ;   ALLOCATE( zwb(jpoce,jpksed) )
97
98
99      ! Initialization of data for mass balance calculation
100      !---------------------------------------------------
101      fromsed(:,:) = 0.
102      tosed  (:,:) = 0. 
103      rloss  (:,:) = 0.
104
105
106      ! Initiate gap
107      !--------------
108      zgap(:,:) = 0.
109      DO js = 1, jpsol
110         DO jk = 1, jpksed
111            DO ji = 1, jpoce
112               zgap(ji,jk) = zgap(ji,jk) + solcp(ji,jk,js)
113            END DO
114         ENDDO
115      ENDDO
116
117      zgap(1:jpoce,1:jpksed) = 1. - zgap(1:jpoce,1:jpksed)   
118
119
120      ! Initiate burial rates
121      !-----------------------
122      zwb(:,:) = 0.
123      DO jk = 2, jpksed
124         zfrac =  dtsed / ( dens * por1(jk) )     
125         DO ji = 1, jpoce
126            zwb(ji,jk) = zfrac * raintg(ji)
127         END DO
128      ENDDO
129
130
131      DO ji = 1, jpoce
132         zwb(ji,2) = zwb(ji,2) - zgap(ji,2) * dz(2)
133      ENDDO
134
135      DO jk = 3, jpksed
136         zfrac = por1(jk-1) / por1(jk)
137         DO ji = 1, jpoce
138            zwb(ji,jk) = zwb(ji,jk-1) * zfrac - zgap(ji,jk) * dz(jk)
139         END DO
140      ENDDO
141
142
143      zrainrf(:,:) = 0.
144      DO ji = 1, jpoce
145         IF( raintg(ji) /= 0. )  &
146            &   zrainrf(ji,:) = rainrg(ji,:) / raintg(ji)
147      ENDDO
148
149
150      ! Computation of full and empty solid fraction in each layer
151      ! for all 'burial' case
152      !----------------------------------------------------------
153
154
155      DO ji = 1, jpoce
156
157         ! computation of total weight fraction in sediment
158         !-------------------------------------------------
159         zfilled(:) = 0.
160         DO js = 1, jpsol
161            DO jk = 2, jpksed
162               zfilled(jk) = zfilled(jk) + solcp(ji,jk,js)
163            ENDDO
164         ENDDO
165         
166         DO js = 1, jpsol
167            DO jk = 2, jpksed
168               zsolcpno(jk,js) = solcp(ji,jk,js) / zfilled(jk)
169            ENDDO
170         ENDDO
171
172         ! burial  3 cases:
173         ! zwb > 0 ==> rain > total rection loss
174         ! zwb = 0 ==> rain = 0
175         ! zwb < 0 ==> rain > 0 and rain < total reaction loss
176         !----------------------------------------------------------------
177
178         IF( zwb(ji,jpksed) > 0. ) THEN
179
180            zfull (jpksed) = zfilled(jpksed)
181            zempty(jpksed) = 1. - zfull(jpksed)
182            DO jk = jpksedm1, 2, -1
183               zfull (jk) = zfilled(jk)
184               zfull (jk) = zfull(jk) - zempty(jk+1) * dvolsp(jk)
185               zempty(jk) = 1. - zfull(jk)
186            ENDDO
187
188            ! Computation of solid sediment species
189            !--------------------------------------
190            ! push entire sediment column downward to account rest of rain
191            DO js = 1, jpsol
192               DO jk = jpksed, 3, -1
193                  solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)
194               ENDDO
195
196               solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)
197
198               DO jk = 2, jpksed
199                  zsolcpno(jk,js) = solcp(ji,jk,js)
200               END DO
201            ENDDO
202
203            zrest = zwb(ji,jpksed) * cpor
204            ! what is remaining is less than dz(2)
205            IF( zrest <= dz(2) ) THEN           
206
207               zfromup(2) = zrest / dz(2)
208               DO jk = 3, jpksed
209                  zfromup(jk) = zwb(ji,jpksed) * ckpor(jk) / dz(jk)
210               ENDDO
211               DO js = 1, jpsol
212                  zfromce = 1. - zfromup(2)
213                  solcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js)
214                  DO jk = 3, jpksed
215                     zfromce = 1. - zfromup(jk)
216                     solcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js)
217                  ENDDO
218                  fromsed(ji,js) = 0.
219                  ! quantities to push in deeper sediment
220                  tosed  (ji,js) = zsolcpno(jpksed,js) &
221                     &           * zwb(ji,jpksed) * dens * por1(jpksed) / mol_wgt(js)
222               ENDDO
223               
224            ELSE ! what is remaining is great than dz(2)
225
226               ntimes = INT( zrest / dz(2) ) + 1
227               IF( ntimes > nztime ) THEN
228                  WRITE( numsed,* ) ' sedadv : rest too large at sediment point ji = ', ji
229                  STOP
230               ENDIF
231               zraipush(1) = dz(2)
232               zrest = zrest - zraipush(1)
233               DO jn = 2, ntimes
234                  IF( zrest >= dz(2) ) THEN
235                     zraipush(jn) = dz(2)
236                     zrest = zrest - zraipush(jn)
237                  ELSE
238                     zraipush(jn) = zrest
239                     zrest = 0.
240                  ENDIF
241               ENDDO
242
243               DO jn = 1, ntimes
244                  DO js = 1, jpsol
245                     DO jk = 2, jpksed
246                        zsolcpno(jk,js) = solcp(ji,jk,js)
247                     END DO
248                  ENDDO
249                 
250                  zfromup(2) = zraipush(jn) / dz(2)
251                  DO jk = 3, jpksed
252                     zfromup(jk) = ( zraipush(jn) / dz(jk) ) * c2por(jk)
253                  ENDDO
254
255                  DO js = 1, jpsol
256                     zfromce = 1. - zfromup(2)
257                     solcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js)
258                     DO jk = 3, jpksed
259                        zfromce = 1. - zfromup(jk)
260                        solcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js)
261                     ENDDO
262                     fromsed(ji,js) = 0.
263                     tosed  (ji,js) = tosed(ji,js) + zsolcpno(jpksed,js) * zraipush(jn) &
264                        &             * dens * por1(2) / mol_wgt(js)
265                  ENDDO
266               ENDDO
267 
268            ENDIF
269
270         ELSE IF( raintg(ji) < eps ) THEN ! rain  = 0
271!! Nadia    rloss(:,:) = rainrm(:,:)   bug ??????           
272
273            rloss(ji,1:jpsol) = rainrm(ji,1:jpsol)
274
275            zfull (2) = zfilled(2)
276            zempty(2) = 1. - zfull(2)
277            DO jk = 3, jpksed
278               zfull (jk) = zfilled(jk)
279               zfull (jk) = zfull (jk) - zempty(jk-1) * dvolsm(jk)
280               zempty(jk) = 1. - zfull(jk)
281            ENDDO
282
283            ! fill boxes with weight fraction from underlying box
284            DO js = 1, jpsol
285               DO jk = 2, jpksedm1
286                  solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)
287               END DO
288               solcp(ji,jpksed,js) = zsolcpno(jpksed,js) * zfull(jpksed)
289               tosed  (ji,js) = 0.
290               fromsed(ji,js) = 0.
291            ENDDO
292            ! for the last layer, one make go up clay
293            solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.
294            !! C. Heinze      fromsed(ji,jsclay) = zempty(jpksed) * 1. * dens * por1(jpksed) / mol_wgt(jsclay)
295            fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay
296
297         ELSE  ! rain > 0 and rain < total reaction loss
298
299
300            DO jk = 2, jpksed
301               zfull (jk) = zfilled(jk)
302               zempty(jk) = 1. - zfull(jk)
303            ENDDO
304
305            ! Determination of indice of layer - ikwneg - where advection is reversed
306            !------------------------------------------------------------------------
307 iflag:     DO jk = 2, jpksed
308               IF( zwb(ji,jk) < 0.  ) THEN
309                  ikwneg = jk
310                  EXIT iflag
311               ENDIF
312            ENDDO iflag
313
314            ! computation of zfull and zempty
315            ! 3 cases : a/ ikwneg=2, b/ikwneg=3...jpksedm1, c/ikwneg=jpksed   
316            !-------------------------------------------------------------     
317            IF( ikwneg == 2 ) THEN ! advection is reversed in the first sediment layer
318
319               zkwnup = rdtsed(ikwneg) * raintg(ji) / dz(ikwneg)
320               zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg)
321               zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1)
322               zempty(ikwneg+1) = 1. - zfull(ikwneg+1)
323               DO jk = ikwneg+2, jpksed
324                  zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk)
325                  zempty(jk) = 1. - zfull(jk)
326               ENDDO
327               DO js = 1, jpsol
328                  solcp(ji,2,js) = zfull(2) * zsolcpno(2,js)+ zkwnlo * zsolcpno(3,js) &
329                     &                                      + zkwnup * zrainrf(ji,js)
330                  DO jk = 3, jpksedm1
331                     solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)
332                  ENDDO
333                  solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)
334                  tosed(ji,js)   = 0.
335                  fromsed(ji,js) = 0.
336               ENDDO
337               solcp(ji,jpksed,jsclay) =  solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.
338               !! C. Heinze  fromsed(ji,jsclay) = zempty(jpksed) * 1. * dens * por1(jpksed) / mol_wgt(jsclay)
339               fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay
340               
341            ELSE IF( ikwneg == jpksed ) THEN
342
343               zkwnup = ABS( zwb(ji,ikwneg-1) ) * dvolsm(ikwneg) / dz(ikwneg)
344               zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg)
345               zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1)
346               zempty(ikwneg-1) = 1. - zfull(ikwneg-1)
347               DO jk = ikwneg-2, 2, -1
348                  zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk)
349                  zempty(jk) = 1. - zfull(jk) 
350               ENDDO
351               DO  js = 1, jpsol
352                  solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)
353               ENDDO
354               DO  js = 1, jpsol             
355                  DO jk = jpksedm1, 3, -1
356                     solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)
357                  ENDDO
358                  solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) &
359                     &                       + zkwnup * zsolcpno(jpksedm1,js)
360                  tosed(ji,js)   = 0.
361                  fromsed(ji,js) = 0.
362               ENDDO
363               solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zkwnlo * 1.
364               ! Heinze  fromsed(ji,jsclay) = zkwnlo * 1. * dens * por1(jpksed) / mol_wgt(jsclay)
365               fromsed(ji,jsclay) = zkwnlo * 1.* por1clay
366            ELSE   ! 2 < ikwneg(ji) <= jpksedm1
367
368               zkwnup = ABS( zwb(ji,ikwneg-1) ) * por1(ikwneg-1) / ( dz(ikwneg) * por1(ikwneg) )
369               zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg)
370
371               IF( ikwneg > 3 ) THEN
372
373                  zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1)
374                  zempty(ikwneg-1) = 1. - zfull(ikwneg-1) 
375                  DO jk = ikwneg-2, 2, -1
376                     zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk)
377                     zempty(jk)    = 1. - zfull(jk) 
378                  ENDDO
379                  DO js = 1, jpsol
380                     solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)
381                  ENDDO
382                  DO js = 1, jpsol
383                     DO jk = ikwneg-1, 3, -1
384                        solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)
385                     ENDDO
386                  ENDDO
387               ELSE ! ikw = 3
388
389
390                  zfull (2) = zfilled(2) - zkwnup * dvolsm(3)
391                  zempty(2) = 1. - zfull(2)
392                  DO js = 1, jpsol
393                     solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)
394                  ENDDO
395               ENDIF
396
397               IF( ikwneg < jpksedm1) THEN
398
399                  zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1)
400                  zempty(ikwneg+1) = 1. - zfull(ikwneg+1) 
401                  DO jk = ikwneg+2, jpksed
402                     zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk)
403                     zempty(jk) = 1. - zfull(jk)
404                  ENDDO
405                  DO js = 1, jpsol
406                     DO jk = ikwneg+1, jpksedm1
407                        solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js) 
408                     ENDDO
409                     solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)
410                  ENDDO
411                  solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.
412               ELSE
413
414                  zfull (jpksed) = zfilled(jpksed) - zkwnlo * dvolsm(jpksed)
415                  zempty(jpksed) = 1. - zfull(jpksed)               
416                  DO js = 1, jpsol
417                     solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)
418                  ENDDO
419                  solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.
420               ENDIF ! jpksedm1
421               
422               ! ikwneg = jpksedm1  ; ikwneg+1 = jpksed ; ikwneg-1 = jpksed - 2
423               DO js = 1, jpsol
424                  solcp(ji,ikwneg,js) =  zfull(ikwneg) * zsolcpno(ikwneg  ,js) &
425                     &                +  zkwnup        * zsolcpno(ikwneg-1,js) &
426                     &                +  zkwnlo        * zsolcpno(ikwneg+1,js)   
427                  tosed  (ji,js)   = 0.
428                  fromsed(ji,js)   = 0.
429               ENDDO
430               ! Heinze  fromsed(ji,jsclay) = zempty * 1. * dens * por1(jpksed) / mol_wgt(jsclay)
431               fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay
432            ENDIF ! ikwneg(ji) = 2
433         ENDIF  ! zwb > 0
434      ENDDO  ! ji = 1, jpoce
435
436      rainrm(:,:) = 0.
437      rainrg(:,:) = 0.
438      raintg(:)   = 0.
439
440
441      DEALLOCATE( zgap  ) ;  DEALLOCATE( zwb )
442
443
444   END SUBROUTINE sed_adv
445#else
446   !!======================================================================
447   !! MODULE sedbtb  :   Dummy module
448   !!======================================================================
449CONTAINS
450   SUBROUTINE sed_adv( kt )         ! Empty routine
451      INTEGER, INTENT(in) :: kt
452      WRITE(*,*) 'sed_adv: You should not have seen this print! error?', kt
453   END SUBROUTINE sed_adv
454
455   !!======================================================================
456
457#endif
458END MODULE sedadv
Note: See TracBrowser for help on using the repository browser.