source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedadv.F90 @ 10345

Last change on this file since 10345 was 10345, checked in by smasson, 23 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10344, see #2133

  • Property svn:keywords set to Id
File size: 16.8 KB
Line 
1MODULE sedadv
2   !!======================================================================
3   !!              ***  MODULE  sedadv  ***
4   !!    Sediment : vertical advection and burial
5   !!=====================================================================
6   !! * Modules used
7   !!----------------------------------------------------------------------
8   !!   sed_adv :
9   !!----------------------------------------------------------------------
10   USE sed     ! sediment global variable
11   USE lib_mpp         ! distribued memory computing library
12
13   IMPLICIT NONE
14   PRIVATE
15
16   PUBLIC sed_adv
17   PUBLIC sed_adv_alloc
18
19   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: dvolsp, dvolsm
20   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: 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, nztime, ikwneg
52     
53      REAL(wp), DIMENSION(jpksed,jpsol) :: zsolcpno
54      REAL(wp), DIMENSION(jpksed)       :: zfilled, zfull, zfromup, zempty
55      REAL(wp), DIMENSION(jpoce,jpksed) :: zgap, zwb
56      REAL(wp), DIMENSION(jpoce,jpsol) :: zrainrf
57      REAL(wp), DIMENSION(:  ), ALLOCATABLE :: zraipush
58
59      REAL(wp) :: zkwnup, zkwnlo, zfrac,  zfromce, zrest, sumtot, zsumtot1
60
61      !------------------------------------------------------------------------
62
63
64      IF( ln_timing )  CALL timing_start('sed_adv')
65!
66      IF( kt == nitsed000 ) THEN
67         IF (lwp) THEN
68            WRITE(numsed,*) ' '
69            WRITE(numsed,*) ' sed_adv : vertical sediment advection  '
70            WRITE(numsed,*) ' '
71         ENDIF
72         por1clay = denssol * por1(jpksed) * dz(jpksed)
73         cpor     = por1(jpksed) / por1(2)
74         DO jk = 2, jpksed
75            c2por(jk) = por1(2)      / por1(jk)
76            ckpor(jk) = por1(jpksed) / por1(jk)
77         ENDDO
78         DO jk = jpksedm1, 2, -1
79            dvolsp(jk) = vols(jk+1) / vols(jk)
80         ENDDO
81         DO jk = 3, jpksed
82           dvolsm(jk) = vols(jk-1) / vols(jk)
83        ENDDO
84      ENDIF
85
86      ! Initialization of data for mass balance calculation
87      !---------------------------------------------------
88      fromsed(:,:) = 0.
89      tosed  (:,:) = 0. 
90      rloss  (:,:) = 0.
91      ikwneg = 1
92      nztime = jpksed
93
94      ALLOCATE( zraipush(nztime) )
95
96      ! Initiate gap
97      !--------------
98      zgap(:,:) = 0.
99      DO js = 1, jpsol
100         DO jk = 1, jpksed
101            DO ji = 1, jpoce
102               zgap(ji,jk) = zgap(ji,jk) + solcp(ji,jk,js)
103            END DO
104         ENDDO
105      ENDDO
106
107      zgap(1:jpoce,1:jpksed) = 1. - zgap(1:jpoce,1:jpksed)   
108
109      ! Initiate burial rates
110      !-----------------------
111      zwb(:,:) = 0.
112      DO jk = 2, jpksed
113         zfrac =  dtsed / ( denssol * por1(jk) )     
114         DO ji = 1, jpoce
115            zwb(ji,jk) = zfrac * raintg(ji)
116         END DO
117      ENDDO
118
119
120      DO ji = 1, jpoce
121         zwb(ji,2) = zwb(ji,2) - zgap(ji,2) * dz(2)
122      ENDDO
123
124      DO jk = 3, jpksed
125         zfrac = por1(jk-1) / por1(jk)
126         DO ji = 1, jpoce
127            zwb(ji,jk) = zwb(ji,jk-1) * zfrac - zgap(ji,jk) * dz(jk)
128         END DO
129      ENDDO
130
131      zrainrf(:,:) = 0.
132      DO ji = 1, jpoce
133         IF( raintg(ji) /= 0. )  &
134            &   zrainrf(ji,:) = rainrg(ji,:) / raintg(ji)
135      ENDDO
136
137
138      ! Computation of full and empty solid fraction in each layer
139      ! for all 'burial' case
140      !----------------------------------------------------------
141
142
143      DO ji = 1, jpoce
144
145         ! computation of total weight fraction in sediment
146         !-------------------------------------------------
147         zfilled(:) = 0.
148         DO js = 1, jpsol
149            DO jk = 2, jpksed
150               zfilled(jk) = zfilled(jk) + solcp(ji,jk,js)
151            ENDDO
152         ENDDO
153         
154         DO js = 1, jpsol
155            DO jk = 2, jpksed
156               zsolcpno(jk,js) = solcp(ji,jk,js) / zfilled(jk)
157            ENDDO
158         ENDDO
159
160         ! burial  3 cases:
161         ! zwb > 0 ==> rain > total rection loss
162         ! zwb = 0 ==> rain = 0
163         ! zwb < 0 ==> rain > 0 and rain < total reaction loss
164         !----------------------------------------------------------------
165
166         IF( zwb(ji,jpksed) > 0. ) THEN
167
168            zfull (jpksed) = zfilled(jpksed)
169            zempty(jpksed) = 1. - zfull(jpksed)
170            DO jk = jpksedm1, 2, -1
171               zfull (jk) = zfilled(jk)
172               zfull (jk) = zfull(jk) - zempty(jk+1) * dvolsp(jk)
173               zempty(jk) = 1. - zfull(jk)
174            ENDDO
175
176            ! Computation of solid sediment species
177            !--------------------------------------
178            ! push entire sediment column downward to account rest of rain
179            DO js = 1, jpsol
180               DO jk = jpksed, 3, -1
181                  solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)
182               ENDDO
183
184               solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)
185
186               DO jk = 2, jpksed
187                  zsolcpno(jk,js) = solcp(ji,jk,js)
188               END DO
189            ENDDO
190
191            zrest = zwb(ji,jpksed) * cpor
192            ! what is remaining is less than dz(2)
193            IF( zrest <= dz(2) ) THEN           
194
195               zfromup(2) = zrest / dz(2)
196               DO jk = 3, jpksed
197                  zfromup(jk) = zwb(ji,jpksed) * ckpor(jk) / dz(jk)
198               ENDDO
199               DO js = 1, jpsol
200                  zfromce = 1. - zfromup(2)
201                  solcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js)
202                  DO jk = 3, jpksed
203                     zfromce = 1. - zfromup(jk)
204                     solcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js)
205                  ENDDO
206                  fromsed(ji,js) = 0.
207                  ! quantities to push in deeper sediment
208                  tosed  (ji,js) = zsolcpno(jpksed,js) &
209                     &           * zwb(ji,jpksed) * denssol * por1(jpksed)
210               ENDDO
211
212            ELSE ! what is remaining is great than dz(2)
213
214               ntimes = INT( zrest / dz(2) ) + 1
215               IF( ntimes > nztime ) CALL ctl_stop( 'STOP', 'sed_adv : rest too large ' )
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                        &             * denssol * por1(2) 
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            fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay
280         ELSE  ! rain > 0 and rain < total reaction loss
281
282
283            DO jk = 2, jpksed
284               zfull (jk) = zfilled(jk)
285               zempty(jk) = 1. - zfull(jk)
286            ENDDO
287
288            ! Determination of indice of layer - ikwneg - where advection is reversed
289            !------------------------------------------------------------------------
290 iflag:     DO jk = 2, jpksed
291               IF( zwb(ji,jk) < 0.  ) THEN
292                  ikwneg = jk
293                  EXIT iflag
294               ENDIF
295            ENDDO iflag
296
297            ! computation of zfull and zempty
298            ! 3 cases : a/ ikwneg=2, b/ikwneg=3...jpksedm1, c/ikwneg=jpksed   
299            !-------------------------------------------------------------     
300            IF( ikwneg == 2 ) THEN ! advection is reversed in the first sediment layer
301
302               zkwnup = rdtsed(ikwneg) * raintg(ji) / dz(ikwneg)
303               zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg)
304               zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1)
305               zempty(ikwneg+1) = 1. - zfull(ikwneg+1)
306               DO jk = ikwneg+2, jpksed
307                  zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk)
308                  zempty(jk) = 1. - zfull(jk)
309               ENDDO
310               DO js = 1, jpsol
311                  solcp(ji,2,js) = zfull(2) * zsolcpno(2,js)+ zkwnlo * zsolcpno(3,js) &
312                     &                                      + zkwnup * zrainrf(ji,js)
313                  DO jk = 3, jpksedm1
314                     solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)
315                  ENDDO
316                  solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)
317                  tosed(ji,js)   = 0.
318                  fromsed(ji,js) = 0.
319               ENDDO
320               solcp(ji,jpksed,jsclay) =  solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.
321               !! C. Heinze  fromsed(ji,jsclay) = zempty(jpksed) * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)
322               fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay
323               
324            ELSE IF( ikwneg == jpksed ) THEN
325
326               zkwnup = ABS( zwb(ji,ikwneg-1) ) * dvolsm(ikwneg) / dz(ikwneg)
327               zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg)
328               zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1)
329               zempty(ikwneg-1) = 1. - zfull(ikwneg-1)
330               DO jk = ikwneg-2, 2, -1
331                  zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk)
332                  zempty(jk) = 1. - zfull(jk) 
333               ENDDO
334               DO  js = 1, jpsol
335                  solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)
336               ENDDO
337               DO  js = 1, jpsol
338                  DO jk = jpksedm1, 3, -1
339                     solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)
340                  ENDDO
341                  solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) &
342                     &                       + zkwnup * zsolcpno(jpksedm1,js)
343                  tosed(ji,js)   = 0.
344                  fromsed(ji,js) = 0.
345               ENDDO
346               solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zkwnlo * 1.
347               ! Heinze  fromsed(ji,jsclay) = zkwnlo * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)
348               fromsed(ji,jsclay) = zkwnlo * 1.* por1clay
349            ELSE   ! 2 < ikwneg(ji) <= jpksedm1
350
351               zkwnup = ABS( zwb(ji,ikwneg-1) ) * por1(ikwneg-1) / ( dz(ikwneg) * por1(ikwneg) )
352               zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg)
353
354               IF( ikwneg > 3 ) THEN
355
356                  zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1)
357                  zempty(ikwneg-1) = 1. - zfull(ikwneg-1) 
358                  DO jk = ikwneg-2, 2, -1
359                     zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk)
360                     zempty(jk)    = 1. - zfull(jk) 
361                  ENDDO
362                  DO js = 1, jpsol
363                     solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)
364                  ENDDO
365                  DO js = 1, jpsol
366                     DO jk = ikwneg-1, 3, -1
367                        solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)
368                     ENDDO
369                  ENDDO
370               ELSE ! ikw = 3
371
372
373                  zfull (2) = zfilled(2) - zkwnup * dvolsm(3)
374                  zempty(2) = 1. - zfull(2)
375                  DO js = 1, jpsol
376                     solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)
377                  ENDDO
378               ENDIF
379
380               IF( ikwneg < jpksedm1) THEN
381
382                  zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1)
383                  zempty(ikwneg+1) = 1. - zfull(ikwneg+1) 
384                  DO jk = ikwneg+2, jpksed
385                     zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk)
386                     zempty(jk) = 1. - zfull(jk)
387                  ENDDO
388                  DO js = 1, jpsol
389                     DO jk = ikwneg+1, jpksedm1
390                        solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js) 
391                     ENDDO
392                     solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)
393                  ENDDO
394                  solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.
395               ELSE
396
397                  zfull (jpksed) = zfilled(jpksed) - zkwnlo * dvolsm(jpksed)
398                  zempty(jpksed) = 1. - zfull(jpksed)               
399                  DO js = 1, jpsol
400                     solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)
401                  ENDDO
402                  solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.
403               ENDIF ! jpksedm1
404               
405               ! ikwneg = jpksedm1  ; ikwneg+1 = jpksed ; ikwneg-1 = jpksed - 2
406               DO js = 1, jpsol
407                  solcp(ji,ikwneg,js) =  zfull(ikwneg) * zsolcpno(ikwneg  ,js) &
408                     &                +  zkwnup        * zsolcpno(ikwneg-1,js) &
409                     &                +  zkwnlo        * zsolcpno(ikwneg+1,js)   
410                  tosed  (ji,js)   = 0.
411                  fromsed(ji,js)   = 0.
412               ENDDO
413               ! Heinze  fromsed(ji,jsclay) = zempty * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)
414               fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay
415
416            ENDIF ! ikwneg(ji) = 2
417         ENDIF  ! zwb > 0
418      ENDDO  ! ji = 1, jpoce
419
420      rainrm(:,:) = 0.
421      rainrg(:,:) = 0.
422      raintg(:)   = 0.
423
424      DEALLOCATE( zraipush )
425
426      IF( ln_timing )  CALL timing_stop('sed_adv')
427
428   END SUBROUTINE sed_adv
429
430
431   INTEGER FUNCTION sed_adv_alloc()
432      !!----------------------------------------------------------------------
433      !!                     ***  ROUTINE p4z_prod_alloc  ***
434      !!----------------------------------------------------------------------
435      ALLOCATE( dvolsp(jpksed), dvolsm(jpksed), c2por(jpksed),         &
436      &         ckpor(jpksed) ,           STAT = sed_adv_alloc )
437      !
438      IF( sed_adv_alloc /= 0 ) CALL ctl_warn('sed_adv_alloc : failed to allocate arrays.')
439      !
440   END FUNCTION sed_adv_alloc
441
442END MODULE sedadv
Note: See TracBrowser for help on using the repository browser.