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 trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/SED – NEMO

source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedadv.F90 @ 4528

Last change on this file since 4528 was 3443, checked in by cetlod, 12 years ago

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