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 NEMO/branches/UKMO/dev_r9888_GO6_mixing/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/UKMO/dev_r9888_GO6_mixing/src/TOP/PISCES/SED/sedadv.F90 @ 9900

Last change on this file since 9900 was 9900, checked in by davestorkey, 6 years ago

UKMO dev_r9888_GO6_mixing branch: clear SVN keywords.

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