- Timestamp:
- 2018-11-16T16:13:30+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/dev_r9950_GO6_mixing/src/TOP/PISCES/SED/sedini.F90
r9950 r10323 4 4 !! Sediment : define sediment variables 5 5 !!===================================================================== 6 #if defined key_sed7 6 8 7 !!---------------------------------------------------------------------- … … 11 10 !! * Modules used 12 11 USE sed ! sediment global variable 13 USE seddta 14 USE sedrst 15 USE sedco3 16 USE sedchem 12 USE sed_oce 17 13 USE sedarr 14 USE sedadv 15 USE trcdmp_sed 16 USE trcdta 18 17 USE iom 18 USE lib_mpp ! distribued memory computing library 19 19 20 20 … … 24 24 !! Module variables 25 25 REAL(wp) :: & 26 sisat = 800. , & !: saturation for si [ mol.l-1] 27 claysat = 0. , & !: saturation for clay [ mol.l-1] 26 sedzmin = 0.3 , & !: Minimum vertical spacing 27 sedhmax = 10.0 , & !: Maximum depth of the sediment 28 sedkth = 5.0 , & !: Default parameters 29 sedacr = 3.0 !: Default parameters 30 31 REAL(wp) :: & 32 porsurf = 0.95 , & !: Porosity at the surface 33 porinf = 0.75 , & !: Porosity at infinite depth 34 rhox = 2.0 !: Vertical length scale of porosity variation 35 36 REAL(wp) :: & 28 37 rcopal = 40. , & !: reactivity for si [l.mol-1.an-1] 29 rcclay = 0. , & !: reactivity for clay [l.mol-1.an-1]30 38 dcoef = 8.e-6 !: diffusion coefficient (*por) [cm**2/s] 31 39 40 REAL(wp), PUBLIC :: & 41 redO2 = 172. , & !: Redfield coef for Oxygen 42 redNo3 = 16. , & !: Redfield coef for Nitrate 43 redPo4 = 1. , & !: Redfield coef for Phosphate 44 redC = 122. , & !: Redfield coef for Carbon 45 redfep = 0.175 , & !: Ratio for iron bound phosphorus 46 rcorgl = 50. , & !: reactivity for POC/O2 [l.mol-1.an-1] 47 rcorgs = 0.5 , & !: reactivity of the semi-labile component 48 rcorgr = 1E-4 , & !: reactivity of the refractory component 49 rcnh4 = 10E6 , & !: reactivity for O2/NH4 [l.mol-1.an-1] 50 rch2s = 1.E5 , & !: reactivity for O2/ODU [l.mol-1.an-1] 51 rcfe2 = 5.E8 , & !: reactivity for O2/Fe2+ [l.mol-1.an-1] 52 rcfeh2s = 1.E4 , & !: Reactivity for FEOH/H2S [l.mol-1.an-1] 53 rcfes = 1.E5 , & !: Reactivity for FE2+/H2S [l.mol-1.an-1] 54 rcfeso = 3.E5 , & !: Reactivity for FES/O2 [l.mol-1.an-1] 55 xksedo2 = 5E-6 , & !: half-sturation constant for oxic remin. 56 xksedno3 = 5E-6 , & !: half-saturation constant for denitrification 57 xksedfeo = 0.6 , & !: half-saturation constant for iron remin 58 xksedso4 = 2E-3 !: half-saturation constant for SO4 remin 59 32 60 REAL(wp) :: & 33 redO2 = 172. , & !: Redfield coef for Oxygene 34 redNo3 = 16. , & !: Redfield coef for Nitrates 35 redPo4 = 1. , & !: Redfield coef for Phosphates 36 redC = 122. , & !: Redfield coef for Carbone 37 redDnit = 97.6 , & !: Redfield coef for denitrification 38 rcorg = 50. , & !: reactivity for POC/O2 [l.mol-1.an-1] 39 o2seuil = 1. , & !: threshold O2 concentration for [mol.l-1] 40 rcorgN = 50. !: reactivity for POC/No3 [l.mol-1.an-1] 41 42 REAL(wp) :: & 43 rccal = 1000. !: reactivity for calcite [l.mol-1.an-1] 44 45 REAL(wp) :: & 46 dbiot = 15. !: coefficient for bioturbation [cm**2.(1000an-1)] 47 48 LOGICAL :: & 49 ln_rst_sed = .TRUE. !: initialisation from a restart file or not 50 61 rccal = 1000., & !: reactivity for calcite [l.mol-1.an-1] 62 rcligc = 1.E-4 !: L/C ratio in POC 63 64 REAL(wp), PUBLIC :: dbiot = 15. , & !: coefficient for bioturbation [cm**2.(n-1)] 65 dbtbzsc = 10.0 , & !: Vertical scale of variation. If no variation, mixed layer in the sed [cm] 66 xirrzsc = 2.0 !: Vertical scale of irrigation variation. 51 67 REAL(wp) :: & 52 68 ryear = 365. * 24. * 3600. !: 1 year converted in second 69 70 REAL(wp), DIMENSION(jpwat), PUBLIC :: diff1 71 DATA diff1/4.59E-6, 1.104E-5, 4.81E-6 , 9.78E-6, 3.58E-6, 4.01E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6 / 72 73 REAL(wp), DIMENSION(jpwat), PUBLIC :: diff2 74 DATA diff2/1.74E-7, 4.47E-7, 2.51E-7, 3.89E-7, 1.77E-7, 2.5E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7 / 75 76 53 77 54 78 !! * Routine accessibility … … 76 100 !! ! 06-07 (C. Ethe) Re-organization 77 101 !!---------------------------------------------------------------------- 78 INTEGER :: ji, jj, ikt 79 #if defined key_sed_off 80 INTEGER :: numblt 81 INTEGER :: nummsh 82 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta 83 #endif 102 INTEGER :: ji, jj, ikt, ierr 84 103 !!---------------------------------------------------------------------- 85 104 … … 90 109 CALL ctl_opn( numsed, 'sediment.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 91 110 92 WRITE(numsed,*) 93 WRITE(numsed,*) ' L S C E - C E A' 94 WRITE(numsed,*) ' SEDIMENT model' 95 WRITE(numsed,*) ' version 2.0 (2007) ' 96 WRITE(numsed,*) 97 WRITE(numsed,*) 98 99 100 WRITE(numsed,*) ' sed_init : Initialization of sediment module ' 101 WRITE(numsed,*) ' ' 102 103 ! ! Allocate LOBSTER arrays 104 IF( sed_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' ) 105 111 IF (lwp) THEN 112 WRITE(numsed,*) 113 WRITE(numsed,*) ' PISCES framework' 114 WRITE(numsed,*) ' SEDIMENT model' 115 WRITE(numsed,*) ' version 3.0 (2018) ' 116 WRITE(numsed,*) 117 WRITE(numsed,*) 118 ENDIF 119 120 IF(lwp) WRITE(numsed,*) ' sed_init : Initialization of sediment module ' 121 IF(lwp) WRITE(numsed,*) ' ' 122 123 ! Read sediment Namelist 124 !------------------------- 125 CALL sed_init_nam 126 127 ! Allocate SEDIMENT arrays 128 ierr = sed_alloc() 129 ierr = ierr + sed_oce_alloc() 130 ierr = ierr + sed_adv_alloc() 131 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' ) 106 132 107 133 ! Determination of sediments number of points and allocate global variables 108 #if defined key_sed_off109 ! Reading Netcdf Pisces file for deepest water layer thickness [m]110 ! This data will be used as mask to merdge space in 1D vector111 !----------------------------------------------------------------112 113 CALL iom_open ( 'mesh_mask' , nummsh )114 CALL iom_open ( 'e3tbot' , numblt )115 116 ! mask sediment points for outputs117 CALL iom_get( nummsh, jpdom_data, 'tmask' , tmask )118 CALL iom_get( nummsh, jpdom_data, 'mbathy', sbathy )119 120 ! longitude/latitude values121 CALL iom_get ( nummsh, jpdom_data,'nav_lon', glamt(:,:) )122 CALL iom_get ( nummsh, jpdom_data,'nav_lat', gphit(:,:) )123 124 ! bottom layer thickness125 ALLOCATE( zdta(jpi,jpj) )126 CALL iom_get ( numblt, jpdom_data, 'E3TBOT', zdta(:,:) )127 epkbot(:,:) = 0.128 DO jj = 1, jpj129 DO ji = 1, jpi130 ikt = MAX( INT( sbathy(ji,jj) ), 1 )131 IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = zdta(ji,jj)132 ENDDO133 ENDDO134 135 CALL iom_close( nummsh )136 CALL iom_close( numblt )137 138 DEALLOCATE( zdta )139 #else140 141 134 epkbot(:,:) = 0. 142 135 DO jj = 1, jpj … … 144 137 ikt = mbkt(ji,jj) 145 138 IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_1d(ikt) 139 gdepbot(ji,jj) = gdepw_0(ji,jj,ikt) 146 140 ENDDO 147 141 ENDDO 148 #endif149 150 142 151 143 ! computation of total number of ocean points 152 144 !-------------------------------------------- 153 jpoce = COUNT( epkbot(:,:) > 0. ) 154 145 sedmask = 0. 146 IF ( COUNT( epkbot(:,:) > 0. ) == 0 ) THEN 147 sedmask = 0. 148 ELSE 149 sedmask = 1. 150 ENDIF 151 jpoce = MAX( COUNT( epkbot(:,:) > 0. ) , 1 ) 155 152 156 153 ! Allocate memory size of global variables … … 159 156 ALLOCATE( rainrm(jpoce,jpsol) ) ; ALLOCATE( rainrg(jpoce,jpsol) ) ; ALLOCATE( raintg(jpoce) ) 160 157 ALLOCATE( dzdep(jpoce) ) ; ALLOCATE( iarroce(jpoce) ) ; ALLOCATE( dzkbot(jpoce) ) 161 ALLOCATE( temp(jpoce) ) ; ALLOCATE( salt(jpoce) ) 158 ALLOCATE( zkbot(jpoce) ) ; ALLOCATE( db(jpoce,jpksed) ) 159 ALLOCATE( temp(jpoce) ) ; ALLOCATE( salt(jpoce) ) 160 ALLOCATE( diff(jpoce,jpksed,jpwat ) ) ; ALLOCATE( irrig(jpoce, jpksed) ) 161 ALLOCATE( wacc(jpoce) ) ; ALLOCATE( fecratio(jpoce) ) 162 162 ALLOCATE( press(jpoce) ) ; ALLOCATE( densSW(jpoce) ) 163 ALLOCATE( hipor(jpoce,jpksed) ) ; ALLOCATE( co3por(jpoce,jpksed) ) 163 ALLOCATE( hipor(jpoce,jpksed) ) ; ALLOCATE( co3por(jpoce,jpksed) ) 164 164 ALLOCATE( dz3d(jpoce,jpksed) ) ; ALLOCATE( volw3d(jpoce,jpksed) ) ; ALLOCATE( vols3d(jpoce,jpksed) ) 165 ALLOCATE( sedligand(jpoce, jpksed) ) 165 166 166 167 ! Initialization of global variables 167 pwcp (:,:,:) = 0. ; pwcp0 (:,:,:) = 0. ; pwcp_dta (:,:) = 0. 168 solcp (:,:,:) = 0. ; solcp0(:,:,:) = 0. ; rainrm_dta(:,:) = 0. 169 rainrm(:,: ) = 0. ; rainrg(:,: ) = 0. ; raintg (: ) = 0. 170 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0. 171 temp (: ) = 0. ; salt (: ) = 0. 172 press (: ) = 0. ; densSW (: ) = 0. 173 hipor (:,: ) = 0. ; co3por (:,: ) = 0. 174 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0. 168 pwcp (:,:,:) = 0. ; pwcp0 (:,:,:) = 0. ; pwcp_dta (:,:) = 0. 169 solcp (:,:,:) = 0. ; solcp0(:,:,:) = 0. ; rainrm_dta(:,:) = 0. 170 rainrm(:,: ) = 0. ; rainrg(:,: ) = 0. ; raintg (: ) = 0. 171 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0. 172 temp (: ) = 0. ; salt (: ) = 0. ; zkbot (: ) = 0. 173 press (: ) = 0. ; densSW (: ) = 0. ; db (:,:) = 0. 174 hipor (:,: ) = 0. ; co3por (:,: ) = 0. ; irrig (:,:) = 0. 175 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0. 176 fecratio(:) = 1E-5 177 sedligand(:,:) = 0.6E-9 175 178 176 179 ! Chemical variables … … 178 181 ALLOCATE( ak1ps (jpoce) ) ; ALLOCATE( ak2ps (jpoce) ) ; ALLOCATE( ak3ps (jpoce) ) ; ALLOCATE( aksis (jpoce) ) 179 182 ALLOCATE( aksps (jpoce) ) ; ALLOCATE( ak12s (jpoce) ) ; ALLOCATE( ak12ps(jpoce) ) ; ALLOCATE( ak123ps(jpoce) ) 180 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) 183 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) ; ALLOCATE( sieqs (jpoce) ) 184 ALLOCATE( aks3s(jpoce) ) ; ALLOCATE( akf3s(jpoce) ) ; ALLOCATE( sulfats(jpoce) ) 185 ALLOCATE( fluorids(jpoce) ) 181 186 182 187 akbs (:) = 0. ; ak1s (:) = 0. ; ak2s (:) = 0. ; akws (:) = 0. 183 188 ak1ps (:) = 0. ; ak2ps (:) = 0. ; ak3ps (:) = 0. ; aksis (:) = 0. 184 189 aksps (:) = 0. ; ak12s (:) = 0. ; ak12ps(:) = 0. ; ak123ps(:) = 0. 185 borats(:) = 0. ; calcon2(:) = 0. 190 borats(:) = 0. ; calcon2(:) = 0. ; sieqs (:) = 0. 191 aks3s(:) = 0. ; akf3s(:) = 0. ; sulfats(:) = 0. ; fluorids(:) = 0. 186 192 187 193 ! Mass balance calculation … … 191 197 fromsed(:,:) = 0. ; tosed(:,:) = 0. ; rloss(:,:) = 0. ; tokbot(:,:) = 0. 192 198 193 ! Read sediment Namelist194 !-------------------------195 CALL sed_init_nam196 197 199 ! Initialization of sediment geometry 198 200 !------------------------------------ 199 201 CALL sed_init_geom 200 202 201 202 ! sets initial sediment composition 203 ! ( only clay or reading restart file ) 204 !--------------------------------------- 205 CALL sed_init_data 206 207 208 CALL sed_init_wri 209 203 ! Offline specific mode 204 ! --------------------- 205 ln_sediment_offline = .FALSE. 206 207 #if defined key_sed_off 208 ln_sediment_offline = .TRUE. 209 IF (lwp) write(numsed,*) 'Sediment module is run in offline mode' 210 IF (lwp) write(numsed,*) 'key_sed_off is activated at compilation time' 211 IF (lwp) write(numsed,*) 'ln_sed_2way is forced to false' 212 IF (lwp) write(numsed,*) '--------------------------------------------' 213 ln_sed_2way = .FALSE. 214 #endif 215 ! Initialisation of tracer damping 216 ! -------------------------------- 217 IF (ln_sediment_offline) THEN 218 CALL trc_dta_ini(jptra) 219 CALL trc_dmp_sed_ini 220 ENDIF 210 221 211 222 END SUBROUTINE sed_init 212 213 223 214 224 SUBROUTINE sed_init_geom … … 227 237 !! * Modules used 228 238 !! * local declarations 229 INTEGER :: & 230 ji, jj, jk 231 232 #if defined key_sed_off 233 REAL(wp) , DIMENSION (jpi,jpj) :: zdta 234 INTEGER :: numpres 235 #endif 239 INTEGER :: ji, jj, jk, jn 240 REAL(wp) :: za0, za1, zt, zw, zsum, zsur, zprof, zprofw 241 REAL(wp) :: ztmp1, ztmp2 236 242 !---------------------------------------------------------- 237 243 238 WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry '239 WRITE(numsed,*) ' '244 IF(lwp) WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry ' 245 IF(lwp) WRITE(numsed,*) ' ' 240 246 241 247 ! Computation of 1D array of sediments points … … 250 256 END DO 251 257 258 IF ( indoce .EQ. 0 ) THEN 259 indoce = 1 260 iarroce(indoce) = 1 261 ENDIF 262 252 263 IF( indoce .NE. jpoce ) THEN 253 WRITE(numsed,*) ' ' 254 WRITE(numsed,*) 'Warning : number of ocean points indoce = ', indoce, & 255 & ' doesn''t match number of point where epkbot>0 jpoce = ', jpoce 256 WRITE(numsed,*) ' ' 257 WRITE(numsed,*) ' ' 258 STOP 264 CALL ctl_stop( 'STOP', 'sed_ini: number of ocean points indoce doesn''t match number of point' ) 259 265 ELSE 260 WRITE(numsed,*) ' ' 261 WRITE(numsed,*) ' total number of ocean points jpoce = ',jpoce 262 WRITE(numsed,*) ' ' 263 ENDIF 264 265 #if defined key_sed_off 266 267 ! Reading Netcdf Pisces file for deepest water layer thickness [m] 268 ! This data will be used as mask to merdge space in 1D vector 269 !---------------------------------------------------------------- 270 CALL iom_open ( 'pressbot' , numpres ) 271 272 ! pressure in bars 273 CALL iom_get ( numpres, jpdom_data,'BATH', zdta(:,:) ) 274 CALL pack_arr ( jpoce, press(1:jpoce), zdta(1:jpi,1:jpj), iarroce(1:jpoce) ) 275 press(1:jpoce) = press(1:jpoce) * 1.025e-1 276 277 CALL iom_close ( numpres ) 278 #endif 279 280 281 ! mask sediment points for outputs 282 DO jk = 1, jpksed 283 tmasksed(:,:,jk) = tmask(:,:,1) 284 ENDDO 266 IF (lwp) WRITE(numsed,*) ' ' 267 IF (lwp) WRITE(numsed,*) ' total number of ocean points jpoce = ',jpoce 268 IF (lwp) WRITE(numsed,*) ' ' 269 ENDIF 285 270 286 271 ! initialization of dzkbot in [cm] … … 288 273 CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 289 274 dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 275 CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 290 276 291 277 ! Geometry and constants … … 293 279 ! (1st layer= diffusive layer = pur water) 294 280 !------------------------------------------ 295 dz(1) = 0.1 296 dz(2) = 0.3 297 dz(3) = 0.3 298 dz(4) = 0.5 299 dz(5) = 0.5 300 dz(6) = 0.5 301 dz(7) = 1. 302 dz(8) = 1. 303 dz(9) = 1. 304 dz(10) = 2.45 305 dz(11) = 2.45 281 za1 = ( sedzmin - sedhmax / FLOAT(jpksed-1) ) & 282 & / ( TANH((1-sedkth)/sedacr) - sedacr/FLOAT(jpksed-1) * ( LOG( COSH( (jpksed - sedkth) / sedacr) ) & 283 & - LOG( COSH( ( 1 - sedkth) / sedacr) ) ) ) 284 za0 = sedzmin - za1 * TANH( (1-sedkth) / sedacr ) 285 zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr ) ) 286 287 profsedw(1) = 0.0 288 profsed(1) = -dz(1) / 2. 289 DO jk = 2, jpksed 290 zw = REAL( jk , wp ) 291 zt = REAL( jk , wp ) - 0.5_wp 292 profsed(jk) = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) ) ) 293 profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) ) ) 294 END DO 295 296 dz(1) = 0.1 297 DO jk = 2, jpksed 298 dz(jk) = profsedw(jk) - profsedw(jk-1) 299 END DO 306 300 307 301 DO jk = 1, jpksed … … 311 305 ENDDO 312 306 313 ! Depth(jk)= depth of middle of each layer314 !----------------------------------------315 profsed(1) = -dz(1)/ 2.316 DO jk = 2, jpksed317 profsed(jk) = profsed(jk-1) + dz(jk-1) / 2. + dz(jk) / 2.318 ENDDO319 320 307 ! Porosity profile [0] 321 308 !--------------------- 322 por(1) = 1. 323 por(2) = 0.95 324 por(3) = 0.9 325 por(4) = 0.85 326 por(5) = 0.81 327 por(6) = 0.75 328 por(7) = 0.75 329 por(8) = 0.75 330 por(9) = 0.75 331 por(10) = 0.75 332 por(11) = 0.75 309 por(1) = 1.0 310 DO jk = 2, jpksed 311 por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk) ) ) 312 END DO 333 313 334 314 ! inverse of Porosity profile … … 353 333 dz3d(1:jpoce,1) = dz(1) 354 334 355 356 335 !--------------------------------------------- 357 336 ! Molecular weight [g/mol] for solid species 358 337 !--------------------------------------------- 359 360 338 361 339 ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16) … … 371 349 & 0.59 * 27. + 10. * 16. 372 350 351 mol_wgt(jsfeo) = 55.0 + 3.0 * ( 16.0 + 1.0) 352 353 mol_wgt(jsfes) = 55.0 + 32.0 354 373 355 ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1) 374 356 ! But den sity of Poc is an Hydrated material (= POC + 30H2O) … … 377 359 mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ & 378 360 & 16. * 14. + 31. ) / 122. 361 mol_wgt(jspos) = mol_wgt(jspoc) 362 mol_wgt(jspor) = mol_wgt(jspoc) 379 363 380 364 ! CaCO3 … … 384 368 ! Density of solid material in sediment [g/cm**3] 385 369 !------------------------------------------------ 386 dens = 2.6370 denssol = 2.6 387 371 388 372 ! Initialization of diffusion coefficient as function of porosity [cm**2/s] 389 373 !-------------------------------------------------------------------- 390 diff(:) = dcoef * por(:) 374 ! DO jn = 1, jpsol 375 ! DO jk = 1, jpksed 376 ! DO ji = 1, jpoce 377 ! diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) ) 378 ! END DO 379 ! END DO 380 ! END DO 381 382 ! Accumulation rate from Burwicz et al. (2011). This is used to 383 ! compute the flux of clays and minerals 384 ! -------------------------------------------------------------- 385 DO ji = 1, jpoce 386 ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 387 ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 ) 388 wacc(ji) = ztmp1 + ztmp2 389 END DO 391 390 392 391 393 392 ! Initialization of time step as function of porosity [cm**2/s] 394 393 !------------------------------------------------------------------ 395 rdtsed(2:jpksed) = dtsed / ( dens * por1(2:jpksed) )396 397 394 END SUBROUTINE sed_init_geom 398 395 … … 408 405 !!---------------------------------------------------------------------- 409 406 410 INTEGER :: & 411 numnamsed = 28 407 INTEGER :: numnamsed_ref = -1 !! Logical units for namelist sediment 408 INTEGER :: numnamsed_cfg = -1 !! Logical units for namelist sediment 409 INTEGER :: ios ! Local integer output status for namelist read 410 CHARACTER(LEN=20) :: clname 412 411 413 412 TYPE PSED … … 422 421 TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d 423 422 424 NAMELIST/nam_time/nfreq 423 NAMELIST/nam_run/nrseddt,ln_sed_2way 424 NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox 425 425 NAMELIST/nam_trased/sedsol, sedwat 426 426 NAMELIST/nam_diased/seddiag3d, seddiag2d 427 NAMELIST/nam_reac/sisat, claysat, rcopal, rcclay, dcoef 428 NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redDnit, & 429 & rcorg, o2seuil, rcorgN 430 NAMELIST/nam_cal/rccal 431 NAMELIST/nam_dc13/pdb, rc13P, rc13Ca 432 NAMELIST/nam_btb/dbiot 433 NAMELIST/nam_rst/ln_rst_sed 434 435 INTEGER :: jn, jn1 427 NAMELIST/nam_inorg/rcopal, dcoef, rccal, ratligc, rcligc 428 NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs, & 429 & rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfes, rcfeso, & 430 & xksedo2, xksedno3, xksedfeo, xksedso4 431 NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, ln_irrig, xirrzsc 432 NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out 433 434 INTEGER :: ji, jn, jn1 436 435 !------------------------------------------------------- 437 436 438 WRITE(numsed,*) ' sed_init_nam : Read namelists '439 WRITE(numsed,*) ' '437 IF(lwp) WRITE(numsed,*) ' sed_init_nam : Read namelists ' 438 IF(lwp) WRITE(numsed,*) ' ' 440 439 441 440 ! ryear = 1 year converted in second 442 441 !------------------------------------ 443 WRITE(numsed,*) ' ' 444 WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear 445 WRITE(numsed,*) ' ' 442 IF (lwp) THEN 443 WRITE(numsed,*) ' ' 444 WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear 445 WRITE(numsed,*) ' ' 446 ENDIF 446 447 447 448 ! Reading namelist.sed variables 448 449 !--------------------------------- 449 CALL ctl_opn( numnamsed, 'namelist.sediment', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 450 451 dtsed = rdt 450 clname = 'namelist_sediment' 451 IF(lwp) WRITE(numsed,*) ' sed_init_nam : read SEDIMENT namelist' 452 IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~' 453 CALL ctl_opn( numnamsed_ref, TRIM( clname )//'_ref', 'OLD' , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 454 CALL ctl_opn( numnamsed_cfg, TRIM( clname )//'_cfg', 'OLD' , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 455 452 456 nitsed000 = nittrc000 453 457 nitsedend = nitend 454 #if ! defined key_sed_off 455 nwrised = nwritetrc 456 #else 457 nwrised = nwrite 458 #endif 459 ! Diffraction/reaction parameters 460 !---------------------------------- 461 READ( numnamsed, nam_time ) 462 WRITE(numsed,*) ' namelist nam_time' 463 464 #if ! defined key_sed_off 465 nfreq = 1 466 #endif 467 468 WRITE(numsed,*) ' sedimentation time step dtsed = ', dtsed 469 WRITE(numsed,*) ' 1st time step for sediment. nitsed000 = ', nitsed000 470 WRITE(numsed,*) ' last time step for sediment. nitsedend = ', nitsedend 471 WRITE(numsed,*) ' frequency of sediment outputs nwrised = ', nwrised 472 WRITE(numsed,*) ' frequency of restoring inputs data nfreq = ', nfreq 473 WRITE(numsed,*) ' ' 474 475 REWIND( numnamsed ) ! read nattrc 476 READ ( numnamsed, nam_trased ) 458 ! Namelist nam_run 459 REWIND( numnamsed_ref ) ! Namelist nam_run in reference namelist : Pisces variables 460 READ ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901) 461 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist', lwp ) 462 463 REWIND( numnamsed_cfg ) ! Namelist nam_run in reference namelist : Pisces variables 464 READ ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902) 465 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist', lwp ) 466 467 IF (lwp) THEN 468 WRITE(numsed,*) ' namelist nam_run' 469 WRITE(numsed,*) ' Nb of iterations for fast species nrseddt = ', nrseddt 470 WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way 471 ENDIF 472 473 REWIND( numnamsed_ref ) ! Namelist nam_geom in reference namelist : Pisces variables 474 READ ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903) 475 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist', lwp ) 476 477 REWIND( numnamsed_cfg ) ! Namelist nam_geom in reference namelist : Pisces variables 478 READ ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904) 479 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist', lwp ) 480 481 IF (lwp) THEN 482 WRITE(numsed,*) ' namelist nam_geom' 483 WRITE(numsed,*) ' Number of vertical layers jpksed = ', jpksed 484 WRITE(numsed,*) ' Minimum vertical spacing sedzmin = ', sedzmin 485 WRITE(numsed,*) ' Maximum depth of the sediment sedhmax = ', sedhmax 486 WRITE(numsed,*) ' Default parameter sedkth = ', sedkth 487 WRITE(numsed,*) ' Default parameter sedacr = ', sedacr 488 WRITE(numsed,*) ' Sediment porosity at the surface porsurf = ', porsurf 489 WRITE(numsed,*) ' Sediment porosity at infinite depth porinf = ', porinf 490 WRITE(numsed,*) ' Length scale of porosity variation rhox = ', rhox 491 ENDIF 492 493 jpksedm1 = jpksed - 1 494 dtsed = r2dttrc 495 496 REWIND( numnamsed_ref ) ! Namelist nam_trased in reference namelist : Pisces variables 497 READ ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905) 498 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist', lwp ) 499 500 REWIND( numnamsed_cfg ) ! Namelist nam_trased in reference namelist : Pisces variables 501 READ ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906) 502 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist', lwp ) 477 503 478 504 DO jn = 1, jpsol … … 489 515 END DO 490 516 491 WRITE(numsed,*) ' namelist nam_trased' 492 WRITE(numsed,*) ' ' 493 DO jn = 1, jptrased 494 WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn)) 495 WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn)) 496 WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn)) 497 WRITE(numsed,*) ' ' 498 END DO 499 WRITE(numsed,*) ' ' 500 517 IF (lwp) THEN 518 WRITE(numsed,*) ' namelist nam_trased' 519 WRITE(numsed,*) ' ' 520 DO jn = 1, jptrased 521 WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn)) 522 WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn)) 523 WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn)) 524 WRITE(numsed,*) ' ' 525 END DO 526 WRITE(numsed,*) ' ' 527 ENDIF 528 529 REWIND( numnamsed_ref ) ! Namelist nam_diased in reference namelist : Pisces variables 530 READ ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907) 531 907 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist', lwp ) 532 533 REWIND( numnamsed_cfg ) ! Namelist nam_diased in reference namelist : Pisces variables 534 READ ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908) 535 908 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist', lwp ) 501 536 502 REWIND( numnamsed )503 READ( numnamsed, nam_diased )504 505 537 DO jn = 1, jpdia3dsed 506 538 seddia3d(jn) = seddiag3d(jn)%snamesed … … 515 547 END DO 516 548 517 WRITE(numsed,*) ' namelist nam_diased' 518 WRITE(numsed,*) ' ' 519 DO jn = 1, jpdia3dsed 520 WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn)) 521 WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn)) 522 WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn)) 523 WRITE(numsed,*) ' ' 524 END DO 525 526 DO jn = 1, jpdia2dsed 527 WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn)) 528 WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn)) 529 WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn)) 530 WRITE(numsed,*) ' ' 531 END DO 532 533 WRITE(numsed,*) ' ' 534 535 536 ! Diffraction/reaction parameters 549 IF (lwp) THEN 550 WRITE(numsed,*) ' namelist nam_diased' 551 WRITE(numsed,*) ' ' 552 DO jn = 1, jpdia3dsed 553 WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn)) 554 WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn)) 555 WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn)) 556 WRITE(numsed,*) ' ' 557 END DO 558 559 DO jn = 1, jpdia2dsed 560 WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn)) 561 WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn)) 562 WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn)) 563 WRITE(numsed,*) ' ' 564 END DO 565 566 WRITE(numsed,*) ' ' 567 ENDIF 568 569 ! Inorganic chemistry parameters 537 570 !---------------------------------- 538 REWIND( numnamsed ) 539 READ( numnamsed, nam_reac ) 540 WRITE(numsed,*) ' namelist nam_reac' 541 WRITE(numsed,*) ' saturation for si sisat = ', sisat 542 WRITE(numsed,*) ' saturation for clay claysat = ', claysat 543 WRITE(numsed,*) ' reactivity for Si rcopal = ', rcopal 544 WRITE(numsed,*) ' reactivity for clay rcclay = ', rcclay 545 WRITE(numsed,*) ' diff. coef for por. dcoef = ', dcoef 546 WRITE(numsed,*) ' ' 547 571 REWIND( numnamsed_ref ) ! Namelist nam_inorg in reference namelist : Pisces variables 572 READ ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909) 573 909 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist', lwp ) 574 575 REWIND( numnamsed_cfg ) ! Namelist nam_inorg in reference namelist : Pisces variables 576 READ ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910) 577 910 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist', lwp ) 578 579 IF (lwp) THEN 580 WRITE(numsed,*) ' namelist nam_inorg' 581 WRITE(numsed,*) ' reactivity for Si rcopal = ', rcopal 582 WRITE(numsed,*) ' diff. coef for por. dcoef = ', dcoef 583 WRITE(numsed,*) ' reactivity for calcite rccal = ', rccal 584 WRITE(numsed,*) ' L/C ratio in POC ratligc = ', ratligc 585 WRITE(numsed,*) ' reactivity for ligands rcligc = ', rcligc 586 WRITE(numsed,*) ' ' 587 ENDIF 548 588 549 589 ! Unity conversion to get saturation conc. psat in [mol.l-1] 550 590 ! and reactivity rc in [l.mol-1.s-1] 551 591 !---------------------------------------------------------- 552 sat_sil = sisat * 1.e-6 553 sat_clay = claysat * 1.e-6 554 555 reac_sil = rcopal / ryear 556 reac_clay = rcclay / ryear 557 592 reac_sil = rcopal / ryear 593 reac_ligc = rcligc / ryear 558 594 559 595 ! Additional parameter linked to POC/O2/No3/Po4 560 596 !---------------------------------------------- 561 REWIND( numnamsed ) 562 READ( numnamsed, nam_poc ) 563 WRITE(numsed,*) ' namelist nam_poc' 564 WRITE(numsed,*) ' Redfield coef for oxy redO2 = ', redO2 565 WRITE(numsed,*) ' Redfield coef for no3 redNo3 = ', redNo3 566 WRITE(numsed,*) ' Redfield coef for po4 redPo4 = ', redPo4 567 WRITE(numsed,*) ' Redfield coef for carbone redC = ', redC 568 WRITE(numsed,*) ' Redfield coef for denitri redDnit = ', redDnit 569 WRITE(numsed,*) ' reactivity for POC/O2 rcorg = ', rcorg 570 WRITE(numsed,*) ' threshold O2 concen. o2seuil = ', o2seuil 571 WRITE(numsed,*) ' reactivity for POC/NO3 rcorgN = ', rcorgN 572 WRITE(numsed,*) ' ' 573 574 575 so2ut = redO2 / redC 576 srno3 = redNo3 / redC 577 spo4r = redPo4 / redC 578 srDnit = redDnit / redC 579 sthro2 = o2seuil * 1.e-6 ! threshold O2 concen. in [mol.l-1] 597 REWIND( numnamsed_ref ) ! Namelist nam_poc in reference namelist : Pisces variables 598 READ ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911) 599 911 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist', lwp ) 600 601 REWIND( numnamsed_cfg ) ! Namelist nam_poc in reference namelist : Pisces variables 602 READ ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912) 603 912 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist', lwp ) 604 605 IF (lwp) THEN 606 WRITE(numsed,*) ' namelist nam_poc' 607 WRITE(numsed,*) ' Redfield coef for oxy redO2 = ', redO2 608 WRITE(numsed,*) ' Redfield coef for no3 redNo3 = ', redNo3 609 WRITE(numsed,*) ' Redfield coef for po4 redPo4 = ', redPo4 610 WRITE(numsed,*) ' Redfield coef for carbon redC = ', redC 611 WRITE(numsed,*) ' Ration for iron bound P redfep = ', redfep 612 WRITE(numsed,*) ' reactivity for labile POC rcorgl = ', rcorgl 613 WRITE(numsed,*) ' reactivity for semi-refract. POC rcorgs = ', rcorgs 614 WRITE(numsed,*) ' reactivity for refractory POC rcorgr = ', rcorgr 615 WRITE(numsed,*) ' reactivity for NH4 rcnh4 = ', rcnh4 616 WRITE(numsed,*) ' reactivity for H2S rch2s = ', rch2s 617 WRITE(numsed,*) ' reactivity for Fe2+ rcfe2 = ', rcfe2 618 WRITE(numsed,*) ' reactivity for FeOH/H2S rcfeh2s = ', rcfeh2s 619 WRITE(numsed,*) ' reactivity for Fe2+/H2S rcfes = ', rcfes 620 WRITE(numsed,*) ' reactivity for FeS/O2 rcfeso = ', rcfeso 621 WRITE(numsed,*) ' Half-sat. cste for oxic remin xksedo2 = ', xksedo2 622 WRITE(numsed,*) ' Half-sat. cste for denit. xksedno3 = ', xksedno3 623 WRITE(numsed,*) ' Half-sat. cste for iron remin xksedfeo = ', xksedfeo 624 WRITE(numsed,*) ' Half-sat. cste for SO4 remin xksedso4 = ', xksedso4 625 WRITE(numsed,*) ' ' 626 ENDIF 627 628 629 so2ut = redO2 / redC 630 srno3 = redNo3 / redC 631 spo4r = redPo4 / redC 632 srDnit = ( (redO2 + 32. ) * 0.8 - redNo3 - redNo3 * 0.6 ) / redC 580 633 ! reactivity rc in [l.mol-1.s-1] 581 reac_poc = rcorg / ryear 582 reac_no3 = rcorgN / ryear 583 584 585 ! Carbonate parameters 586 !--------------------- 587 READ( numnamsed, nam_cal ) 588 WRITE(numsed,*) ' namelist nam_cal' 589 WRITE(numsed,*) ' reactivity for calcite rccal = ', rccal 590 WRITE(numsed,*) ' ' 634 reac_pocl = rcorgl / ryear 635 reac_pocs = rcorgs / ryear 636 reac_pocr = rcorgr / ryear 637 reac_nh4 = rcnh4 / ryear 638 reac_h2s = rch2s / ryear 639 reac_fe2 = rcfe2 / ryear 640 reac_feh2s = rcfeh2s/ ryear 641 reac_fes = rcfes / ryear 642 reac_feso = rcfeso / ryear 591 643 592 644 ! reactivity rc in [l.mol-1.s-1] 593 645 reac_cal = rccal / ryear 594 646 595 596 ! C13 parameters597 !----------------598 READ( numnamsed, nam_dc13 )599 WRITE(numsed,*) ' namelist nam_dc13 '600 WRITE(numsed,*) ' 13C/12C in PD Belemnite PDB = ', pdb601 WRITE(numsed,*) ' 13C/12C in POC = rc13P*PDB rc13P = ', rc13P602 WRITE(numsed,*) ' 13C/12C in CaCO3 = rc13Ca*PDB rc13Ca = ', rc13Ca603 WRITE(numsed,*) ' '604 605 606 647 ! Bioturbation parameter 607 648 !------------------------ 608 READ( numnamsed, nam_btb ) 609 WRITE(numsed,*) ' namelist nam_btb ' 610 WRITE(numsed,*) ' coefficient for bioturbation dbiot = ', dbiot 611 WRITE(numsed,*) ' ' 612 613 ! Unity convertion to get bioturb coefficient in [cm2.s-1] 614 db = dbiot / ( ryear * 1000. ) 649 REWIND( numnamsed_ref ) ! Namelist nam_btb in reference namelist : Pisces variables 650 READ ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913) 651 913 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist', lwp ) 652 653 REWIND( numnamsed_cfg ) ! Namelist nam_btb in reference namelist : Pisces variables 654 READ ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914) 655 914 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist', lwp ) 656 657 IF (lwp) THEN 658 WRITE(numsed,*) ' namelist nam_btb ' 659 WRITE(numsed,*) ' coefficient for bioturbation dbiot = ', dbiot 660 WRITE(numsed,*) ' Depth varying bioturbation ln_btbz = ', ln_btbz 661 WRITE(numsed,*) ' coefficient for btb attenuation dbtbzsc = ', dbtbzsc 662 WRITE(numsed,*) ' Adsorption coefficient of NH4 adsnh4 = ', adsnh4 663 WRITE(numsed,*) ' Bioirrigation in sediment ln_irrig = ', ln_irrig 664 WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc = ', xirrzsc 665 WRITE(numsed,*) ' ' 666 ENDIF 615 667 616 668 ! Initial value (t=0) for sediment pore water and solid components 617 669 !---------------------------------------------------------------- 618 READ( numnamsed, nam_rst ) 619 WRITE(numsed,*) ' namelist nam_rst ' 620 WRITE(numsed,*) ' boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 621 WRITE(numsed,*) ' ' 622 623 CLOSE( numnamsed ) 670 REWIND( numnamsed_ref ) ! Namelist nam_rst in reference namelist : Pisces variables 671 READ ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915) 672 915 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist', lwp ) 673 674 REWIND( numnamsed_cfg ) ! Namelist nam_rst in reference namelist : Pisces variables 675 READ ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916) 676 916 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist', lwp ) 677 678 IF (lwp) THEN 679 WRITE(numsed,*) ' namelist nam_rst ' 680 WRITE(numsed,*) ' boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 681 WRITE(numsed,*) ' ' 682 ENDIF 683 nn_dtsed = nn_dttrc 684 685 CLOSE( numnamsed_cfg ) 686 CLOSE( numnamsed_ref ) 624 687 625 688 END SUBROUTINE sed_init_nam 626 689 627 SUBROUTINE sed_init_data628 !!----------------------------------------------------------------------629 !! *** ROUTINE sed_init_data ***630 !!631 !! ** Purpose : Initialization of sediment module632 !! - sets initial sediment composition633 !! ( only clay or reading restart file )634 !!635 !! History :636 !! ! 06-07 (C. Ethe) original637 !!----------------------------------------------------------------------638 639 ! local variables640 INTEGER :: &641 ji, jk, zhipor642 643 !--------------------------------------------------------------------644 645 646 IF( .NOT. ln_rst_sed ) THEN647 648 WRITE(numsed,*) ' Initilization of default values of sediment components'649 650 ! default values for initial pore water concentrations [mol/l]651 pwcp(:,:,:) = 0.652 ! default value for initial solid component (fraction of dry weight dim=[0])653 ! clay654 solcp(:,:,:) = 0.655 solcp(:,2:jpksed,jsclay) = 1.656 657 ! Initialization of [h+] and [co3--]658 659 zhipor = 8.660 ! Initialization of [h+] in mol/kg661 DO jk = 1, jpksed662 DO ji = 1, jpoce663 hipor (ji,jk) = 10.**( -1. * zhipor )664 ENDDO665 ENDDO666 667 co3por(:,:) = 0.668 669 ELSE670 671 WRITE(numsed,*) ' Initilization of Sediment components from restart'672 673 CALL sed_rst_read674 675 ENDIF676 677 678 ! Load initial Pisces Data for bot. wat. Chem and fluxes679 CALL sed_dta ( nitsed000 )680 681 ! Initialization of chemical constants682 CALL sed_chem ( nitsed000 )683 684 ! Stores initial sediment data for mass balance calculation685 pwcp0 (1:jpoce,1:jpksed,1:jpwat ) = pwcp (1:jpoce,1:jpksed,1:jpwat )686 solcp0(1:jpoce,1:jpksed,1:jpsol ) = solcp(1:jpoce,1:jpksed,1:jpsol)687 688 ! Conversion of [h+] in mol/Kg to get it in mol/l ( multiplication by density)689 DO jk = 1, jpksed690 hipor(1:jpoce,jk) = hipor(1:jpoce,jk) * densSW(1:jpoce)691 ENDDO692 693 694 ! In default case - no restart - sedco3 is run to initiate [h+] and [co32-]695 ! Otherwise initiate values of pH and co3 read in restart696 IF( .NOT. ln_rst_sed ) THEN697 ! sedco3 is run to initiate[h+] [co32-] in mol/l of solution698 CALL sed_co3 ( nitsed000 )699 700 ENDIF701 702 END SUBROUTINE sed_init_data703 704 SUBROUTINE sed_init_wri705 706 INTEGER :: jk707 708 WRITE(numsed,*)' '709 WRITE(numsed,*)'======== Write summary of initial state ============'710 WRITE(numsed,*)' '711 WRITE(numsed,*)' '712 WRITE(numsed,*)'-------------------------------------------------------------------'713 WRITE(numsed,*)' Initial Conditions '714 WRITE(numsed,*)'-------------------------------------------------------------------'715 WRITE(numsed,*)'dzm = dzkbot minimum to calculate ', 0.716 WRITE(numsed,*)'Local zone : jpi, jpj : ',jpi, jpj717 WRITE(numsed,*)'jpoce = ',jpoce,' nbtot pts = ',jpij,' nb earth pts = ',jpij - jpoce718 WRITE(numsed,*)'sublayer thickness dz(1) [cm] : ', dz(1)719 WRITE(numsed,*)'Coeff diff for k=1 (cm2/s) : ',diff(1)720 WRITE(numsed,*)' nb solid comp : ',jpsol721 WRITE(numsed,*)'(1=opal,2=clay,3=POC,4=CaCO3)'722 WRITE(numsed,*)'weight mol 1,2,3,4'723 WRITE(numsed,'(4(F0.2,3X))')mol_wgt(jsopal),mol_wgt(jsclay),mol_wgt(jspoc),mol_wgt(jscal)724 WRITE(numsed,*)'nb dissolved comp',jpwat725 WRITE(numsed,*)'(1=silicic acid,2="dissolved" clay,3=O2,4=DIC,5=Nitrate,&726 &6=Phosphates,7=Alk))'727 WRITE(numsed,*)'Psat (umol/l) for silicic Acid and "dissolved" clay'728 WRITE(numsed,'(2(F0.2,3X))') sat_sil * 1e+6, sat_clay * 1e+6729 WRITE(numsed,*)'reaction rate rc for Op/si,Clay,POC/O2,caco3, POC/No3 (an-1)'730 WRITE(numsed,'(5(F0.2,3X))') reac_sil * ryear, reac_clay * ryear, reac_poc * ryear, &731 reac_cal * ryear, reac_no3 * ryear732 WRITE(numsed,*)'redfield coef C,O,N P Dit '733 WRITE(numsed,'(5(F0.2,3X))')1./spo4r,so2ut/spo4r,srno3/spo4r,spo4r/spo4r,srDnit/spo4r734 WRITE(numsed,*)'threshold for stating denitrification [mol/l]'735 WRITE(numsed,'(1PE8.2)') sthrO2736 WRITE(numsed,*)'-------------------------------------------------------------------'737 WRITE(numsed,*)'Min-Max-Mean'738 WRITE(numsed,*)'For each variable : min, max, moy value observed on selected local zone'739 WRITE(numsed,*)'-------------------------------------------------------------------'740 WRITE(numsed,*)'thickness of the last wet layer dzkbot [m]'741 WRITE(numsed,'(3(F0.2,3X))') MINVAL(dzkbot(1:jpoce))/100.,MAXVAL(dzkbot(1:jpoce))/100.,&742 &SUM(dzkbot(1:jpoce))/jpoce/100.743 WRITE(numsed,*)'temp [°C]'744 WRITE(numsed,'(3(F0.2,3X))') MINVAL(temp(1:jpoce)),MAXVAL(temp(1:jpoce)),&745 & SUM(temp(1:jpoce))/jpoce746 WRITE(numsed,*)'salt o/oo'747 WRITE(numsed,'(3(F0.2,3X))')MINVAL(salt(1:jpoce)),MAXVAL(salt(1:jpoce)),&748 & SUM(salt(1:jpoce))/jpoce749 #if defined key_sed_off750 WRITE(numsed,*)'pressure [bar] (depth in m is about 10*pressure)'751 WRITE(numsed,'(3(F0.2,3X))') MINVAL(press(1:jpoce)),MAXVAL(press(1:jpoce)),&752 & SUM(press(1:jpoce))/jpoce753 #endif754 WRITE(numsed,*)'density of Sea Water'755 WRITE(numsed,'(3(F0.2,3X))') MINVAL(densSW(1:jpoce)), MAXVAL(densSW(1:jpoce)),&756 & SUM(densSW(1:jpoce))/jpoce757 WRITE(numsed,*)''758 WRITE(numsed,*)' Dissolved Components '759 WRITE(numsed,*)' ====================='760 WRITE(numsed,*)'[Si(OH)4] dissolved (1)(k=1)(µmol/l)(and min value in mol/kg of solution)'761 WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwsil))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwsil))*1.e+6,&762 & SUM(pwcp(1:jpoce,1,jwsil))*1.e+6/jpoce,&763 & MINVAL(pwcp(1:jpoce,1,jwsil)*1.e+6/densSW(1:jpoce))764 WRITE(numsed,*)'[O2] dissolved (3) (k=1)(µmol/l)(and min value in mol/kg of solution)'765 WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwoxy))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwoxy))*1.e+6,&766 &SUM(pwcp(1:jpoce,1,jwoxy))*1.e+6/jpoce,&767 &MINVAL(pwcp(1:jpoce,1,jwoxy)*1.e+6/densSW(1:jpoce))768 WRITE(numsed,*)'[DIC] dissolved (4) (k=1)(µmol/l)(and min value in mol/kg of solution)'769 WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwdic))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwdic))*1.e+6,&770 &SUM(pwcp(1:jpoce,1,jwdic))*1.e+6/jpoce,&771 &MINVAL(pwcp(1:jpoce,1,jwdic)*1.e+6/densSW(1:jpoce))772 WRITE(numsed,*)'[NO3] dissolved (5) (k=1)(µmol/l)(and min value in mol/kg of solution)'773 WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwno3))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwno3))*1.e+6,&774 &SUM(pwcp(1:jpoce,1,jwno3))*1.e+6/jpoce,&775 &MINVAL(pwcp(1:jpoce,1,jwno3)*1.e+6/densSW(1:jpoce))776 WRITE(numsed,*)'[PO4] dissolved (6) (k=1)(µmol/l)(and min value in mol/kg of solution)'777 WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwpo4))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwpo4))*1.e+6,&778 &SUM(pwcp(1:jpoce,1,jwpo4))*1.e+6/jpoce,&779 &MINVAL(pwcp(1:jpoce,1,jwpo4)*1.e+6/densSW(1:jpoce))780 WRITE(numsed,*)'[Alk] dissolved (7) (k=1)(µequi)(and min value in mol/kg of solution)'781 WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwalk))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwalk))*1.e+6,&782 &SUM(pwcp(1:jpoce,1,jwalk))*1.e+6/jpoce,&783 &MINVAL(pwcp(1:jpoce,1,jwalk)*1.e+6/densSW(1:jpoce))784 WRITE(numsed,*)'[DIC13] dissolved (8) (k=1)(µmol/l)(and min value in mol/kg of solution)'785 WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwc13))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwc13))*1.e+6,&786 &SUM(pwcp(1:jpoce,1,jwc13))*1.e+6/jpoce,&787 &MINVAL(pwcp(1:jpoce,1,jwc13)*1.e+6/densSW(1:jpoce))788 WRITE(numsed,*)''789 WRITE(numsed,*)' Solid Components '790 WRITE(numsed,*)' ====================='791 WRITE(numsed,*)'nmole of Opale rained per dt'792 WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsopal))*dtsed,MAXVAL(rainrm(1:jpoce,jsopal))*dtsed,&793 &SUM(rainrm(1:jpoce,1))*dtsed/jpoce794 WRITE(numsed,*)'nmole of Clay rained per dt'795 WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsclay))*dtsed,MAXVAL(rainrm(1:jpoce,jsclay))*dtsed,&796 &SUM(rainrm(1:jpoce,jsclay))*dtsed/jpoce797 WRITE(numsed,*)'nmole of POC rained per dt'798 WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jspoc))*dtsed,MAXVAL(rainrm(1:jpoce,jspoc))*dtsed,&799 &SUM(rainrm(1:jpoce,jspoc))*dtsed/jpoce800 WRITE(numsed,*)'nmole of CaCO3 rained per dt'801 WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jscal))*dtsed,MAXVAL(rainrm(1:jpoce,jscal))*dtsed,&802 &SUM(rainrm(1:jpoce,jscal))*dtsed/jpoce803 WRITE(numsed,*)' '804 WRITE(numsed,*)'Weight frac of opal rained (%) '805 WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,&806 &MAXVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,&807 & SUM(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100./jpoce808 WRITE(numsed,*)'Weight frac of clay rained (%) '809 WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,&810 &MAXVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,&811 &SUM(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100./jpoce812 WRITE(numsed,*)'Weight frac of POC rained (%)'813 WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,&814 &MAXVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,&815 &SUM(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100./jpoce816 WRITE(numsed,*)'Weight frac of CaCO3 rained (%)'817 WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,&818 &MAXVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,&819 &SUM(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100./jpoce820 WRITE(numsed,*)''821 WRITE(numsed,*)' Thickness of sediment layer added by particule rain, dzdep cm '822 WRITE(numsed,*)' =============================================================='823 WRITE(numsed,'(3(F0.5,2X))') MINVAL(dzdep(1:jpoce)),MAXVAL(dzdep(1:jpoce)),SUM(dzdep(1:jpoce))/jpoce824 WRITE(numsed,*)''825 WRITE(numsed,*)' chemical constants K1,pK1,K2,pK2,Kw,pKw and Kb pKb (min max) [mol/kgsol] or [mol/kgsol]2 '826 WRITE(numsed,*)' ========================================================================================='827 WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak1s(1:jpoce)),MAXVAL(ak1s(1:jpoce)),-LOG10(MINVAL(ak1s(1:jpoce))),&828 &-LOG10(MAXVAL(ak1s(1:jpoce)))829 WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak2s(1:jpoce)),MAXVAL(ak2s(1:jpoce)),-LOG10(MINVAL(ak2s(1:jpoce))),&830 &-LOG10(MAXVAL(ak2s(1:jpoce)))831 WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akws(1:jpoce)),MAXVAL(akws(1:jpoce)),-LOG10(MINVAL(akws(1:jpoce))),&832 &-LOG10(MAXVAL(akws(1:jpoce)))833 WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akbs(1:jpoce)),MAXVAL(akbs(1:jpoce)),-LOG10(MINVAL(akbs(1:jpoce))),&834 &-LOG10(MAXVAL(akbs(1:jpoce)))835 WRITE(numsed,*)'and boron'836 WRITE(numsed,'(2(1PE10.3,2X))')MINVAL(borats(1:jpoce)),MAXVAL(borats(1:jpoce))837 WRITE(numsed,*)''838 WRITE(numsed,*)' Compo of initial sediment for point jpoce'839 WRITE(numsed,*)' ========================================='840 WRITE(numsed,*)'solcp(1), solcp(2), solcp(3), solcp(4), hipor, pH, co3por'841 DO jk = 1,jpksed842 WRITE(numsed,'(7(1PE10.3,2X))')solcp(jpoce,jk,jsopal),solcp(jpoce,jk,jsclay),solcp(jpoce,jk,jspoc),solcp(jpoce,jk,jscal),&843 & hipor(jpoce,jk),-LOG10(hipor(jpoce,jk)/densSW(jpoce)),co3por(jpoce,jk)844 ENDDO845 WRITE(numsed,'(A82)')'pwcp(1), pwcp(2), pwcp(3), pwcp(4), pwcp(5), pwcp(6), pwcp(7)'846 DO jk = 1, jpksed847 WRITE(numsed,'(7(1PE10.3,2X))')pwcp(jpoce,jk,jwsil),pwcp(jpoce,jk,jwoxy),pwcp(jpoce,jk,jwdic),&848 & pwcp(jpoce,jk,jwno3),pwcp(jpoce,jk,jwpo4),pwcp(jpoce,jk,jwalk),pwcp(jpoce,jk,jwc13)849 ENDDO850 WRITE(numsed,*) ' '851 WRITE(numsed,*) ' End Of Initialization '852 WRITE(numsed,*) ' '853 !854 END SUBROUTINE sed_init_wri855 #else856 !!----------------------------------------------------------------------857 !! Dummy module : NO Sediment model858 !!----------------------------------------------------------------------859 !! $Id$860 CONTAINS861 SUBROUTINE sed_ini ! Empty routine862 END SUBROUTINE sed_ini863 #endif864 865 866 690 END MODULE sedini
Note: See TracChangeset
for help on using the changeset viewer.
