source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedadv.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

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