Changeset 10345 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedini.F90
- Timestamp:
- 2018-11-21T11:25:53+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedini.F90
r5215 r10345 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 trc_oce, ONLY : nn_dttrc 16 USE trcdmp_sed 17 USE trcdta 18 18 USE iom 19 USE lib_mpp ! distribued memory computing library 19 20 20 21 … … 24 25 !! Module variables 25 26 REAL(wp) :: & 26 sisat = 800. , & !: saturation for si [ mol.l-1] 27 claysat = 0. , & !: saturation for clay [ mol.l-1] 27 sedzmin = 0.3 , & !: Minimum vertical spacing 28 sedhmax = 10.0 , & !: Maximum depth of the sediment 29 sedkth = 5.0 , & !: Default parameters 30 sedacr = 3.0 !: Default parameters 31 32 REAL(wp) :: & 33 porsurf = 0.95 , & !: Porosity at the surface 34 porinf = 0.75 , & !: Porosity at infinite depth 35 rhox = 2.0 !: Vertical length scale of porosity variation 36 37 REAL(wp) :: & 28 38 rcopal = 40. , & !: reactivity for si [l.mol-1.an-1] 29 rcclay = 0. , & !: reactivity for clay [l.mol-1.an-1]30 39 dcoef = 8.e-6 !: diffusion coefficient (*por) [cm**2/s] 31 40 41 REAL(wp), PUBLIC :: & 42 redO2 = 172. , & !: Redfield coef for Oxygen 43 redNo3 = 16. , & !: Redfield coef for Nitrate 44 redPo4 = 1. , & !: Redfield coef for Phosphate 45 redC = 122. , & !: Redfield coef for Carbon 46 redfep = 0.175 , & !: Ratio for iron bound phosphorus 47 rcorgl = 50. , & !: reactivity for POC/O2 [l.mol-1.an-1] 48 rcorgs = 0.5 , & !: reactivity of the semi-labile component 49 rcorgr = 1E-4 , & !: reactivity of the refractory component 50 rcnh4 = 10E6 , & !: reactivity for O2/NH4 [l.mol-1.an-1] 51 rch2s = 1.E5 , & !: reactivity for O2/ODU [l.mol-1.an-1] 52 rcfe2 = 5.E8 , & !: reactivity for O2/Fe2+ [l.mol-1.an-1] 53 rcfeh2s = 1.E4 , & !: Reactivity for FEOH/H2S [l.mol-1.an-1] 54 rcfes = 1.E5 , & !: Reactivity for FE2+/H2S [l.mol-1.an-1] 55 rcfeso = 3.E5 , & !: Reactivity for FES/O2 [l.mol-1.an-1] 56 xksedo2 = 5E-6 , & !: half-sturation constant for oxic remin. 57 xksedno3 = 5E-6 , & !: half-saturation constant for denitrification 58 xksedfeo = 0.6 , & !: half-saturation constant for iron remin 59 xksedso4 = 2E-3 !: half-saturation constant for SO4 remin 60 32 61 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 62 rccal = 1000., & !: reactivity for calcite [l.mol-1.an-1] 63 rcligc = 1.E-4 !: L/C ratio in POC 64 65 REAL(wp), PUBLIC :: dbiot = 15. , & !: coefficient for bioturbation [cm**2.(n-1)] 66 dbtbzsc = 10.0 , & !: Vertical scale of variation. If no variation, mixed layer in the sed [cm] 67 xirrzsc = 2.0 !: Vertical scale of irrigation variation. 51 68 REAL(wp) :: & 52 69 ryear = 365. * 24. * 3600. !: 1 year converted in second 70 71 REAL(wp), DIMENSION(jpwat), PUBLIC :: diff1 72 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 / 73 74 REAL(wp), DIMENSION(jpwat), PUBLIC :: diff2 75 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 / 76 77 53 78 54 79 !! * Routine accessibility … … 76 101 !! ! 06-07 (C. Ethe) Re-organization 77 102 !!---------------------------------------------------------------------- 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 103 INTEGER :: ji, jj, ikt, ierr 84 104 !!---------------------------------------------------------------------- 85 105 … … 90 110 CALL ctl_opn( numsed, 'sediment.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 91 111 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 112 IF (lwp) THEN 113 WRITE(numsed,*) 114 WRITE(numsed,*) ' PISCES framework' 115 WRITE(numsed,*) ' SEDIMENT model' 116 WRITE(numsed,*) ' version 3.0 (2018) ' 117 WRITE(numsed,*) 118 WRITE(numsed,*) 119 ENDIF 120 121 IF(lwp) WRITE(numsed,*) ' sed_init : Initialization of sediment module ' 122 IF(lwp) WRITE(numsed,*) ' ' 123 124 ! Read sediment Namelist 125 !------------------------- 126 CALL sed_init_nam 127 128 ! Allocate SEDIMENT arrays 129 ierr = sed_alloc() 130 ierr = ierr + sed_oce_alloc() 131 ierr = ierr + sed_adv_alloc() 132 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' ) 106 133 107 134 ! 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 135 epkbot(:,:) = 0. 142 136 DO jj = 1, jpj … … 144 138 ikt = mbkt(ji,jj) 145 139 IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_1d(ikt) 140 gdepbot(ji,jj) = gdepw_0(ji,jj,ikt) 146 141 ENDDO 147 142 ENDDO 148 #endif149 150 143 151 144 ! computation of total number of ocean points 152 145 !-------------------------------------------- 153 jpoce = COUNT( epkbot(:,:) > 0. ) 154 146 sedmask = 0. 147 IF ( COUNT( epkbot(:,:) > 0. ) == 0 ) THEN 148 sedmask = 0. 149 ELSE 150 sedmask = 1. 151 ENDIF 152 jpoce = MAX( COUNT( epkbot(:,:) > 0. ) , 1 ) 155 153 156 154 ! Allocate memory size of global variables … … 159 157 ALLOCATE( rainrm(jpoce,jpsol) ) ; ALLOCATE( rainrg(jpoce,jpsol) ) ; ALLOCATE( raintg(jpoce) ) 160 158 ALLOCATE( dzdep(jpoce) ) ; ALLOCATE( iarroce(jpoce) ) ; ALLOCATE( dzkbot(jpoce) ) 161 ALLOCATE( temp(jpoce) ) ; ALLOCATE( salt(jpoce) ) 159 ALLOCATE( zkbot(jpoce) ) ; ALLOCATE( db(jpoce,jpksed) ) 160 ALLOCATE( temp(jpoce) ) ; ALLOCATE( salt(jpoce) ) 161 ALLOCATE( diff(jpoce,jpksed,jpwat ) ) ; ALLOCATE( irrig(jpoce, jpksed) ) 162 ALLOCATE( wacc(jpoce) ) ; ALLOCATE( fecratio(jpoce) ) 162 163 ALLOCATE( press(jpoce) ) ; ALLOCATE( densSW(jpoce) ) 163 ALLOCATE( hipor(jpoce,jpksed) ) ; ALLOCATE( co3por(jpoce,jpksed) ) 164 ALLOCATE( hipor(jpoce,jpksed) ) ; ALLOCATE( co3por(jpoce,jpksed) ) 164 165 ALLOCATE( dz3d(jpoce,jpksed) ) ; ALLOCATE( volw3d(jpoce,jpksed) ) ; ALLOCATE( vols3d(jpoce,jpksed) ) 166 ALLOCATE( sedligand(jpoce, jpksed) ) 165 167 166 168 ! 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. 169 pwcp (:,:,:) = 0. ; pwcp0 (:,:,:) = 0. ; pwcp_dta (:,:) = 0. 170 solcp (:,:,:) = 0. ; solcp0(:,:,:) = 0. ; rainrm_dta(:,:) = 0. 171 rainrm(:,: ) = 0. ; rainrg(:,: ) = 0. ; raintg (: ) = 0. 172 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0. 173 temp (: ) = 0. ; salt (: ) = 0. ; zkbot (: ) = 0. 174 press (: ) = 0. ; densSW (: ) = 0. ; db (:,:) = 0. 175 hipor (:,: ) = 0. ; co3por (:,: ) = 0. ; irrig (:,:) = 0. 176 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0. 177 fecratio(:) = 1E-5 178 sedligand(:,:) = 0.6E-9 175 179 176 180 ! Chemical variables … … 178 182 ALLOCATE( ak1ps (jpoce) ) ; ALLOCATE( ak2ps (jpoce) ) ; ALLOCATE( ak3ps (jpoce) ) ; ALLOCATE( aksis (jpoce) ) 179 183 ALLOCATE( aksps (jpoce) ) ; ALLOCATE( ak12s (jpoce) ) ; ALLOCATE( ak12ps(jpoce) ) ; ALLOCATE( ak123ps(jpoce) ) 180 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) 184 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) ; ALLOCATE( sieqs (jpoce) ) 185 ALLOCATE( aks3s(jpoce) ) ; ALLOCATE( akf3s(jpoce) ) ; ALLOCATE( sulfats(jpoce) ) 186 ALLOCATE( fluorids(jpoce) ) 181 187 182 188 akbs (:) = 0. ; ak1s (:) = 0. ; ak2s (:) = 0. ; akws (:) = 0. 183 189 ak1ps (:) = 0. ; ak2ps (:) = 0. ; ak3ps (:) = 0. ; aksis (:) = 0. 184 190 aksps (:) = 0. ; ak12s (:) = 0. ; ak12ps(:) = 0. ; ak123ps(:) = 0. 185 borats(:) = 0. ; calcon2(:) = 0. 191 borats(:) = 0. ; calcon2(:) = 0. ; sieqs (:) = 0. 192 aks3s(:) = 0. ; akf3s(:) = 0. ; sulfats(:) = 0. ; fluorids(:) = 0. 186 193 187 194 ! Mass balance calculation … … 191 198 fromsed(:,:) = 0. ; tosed(:,:) = 0. ; rloss(:,:) = 0. ; tokbot(:,:) = 0. 192 199 193 ! Read sediment Namelist194 !-------------------------195 CALL sed_init_nam196 197 200 ! Initialization of sediment geometry 198 201 !------------------------------------ 199 202 CALL sed_init_geom 200 203 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 204 ! Offline specific mode 205 ! --------------------- 206 ln_sediment_offline = .FALSE. 207 208 #if defined key_sed_off 209 ln_sediment_offline = .TRUE. 210 IF (lwp) write(numsed,*) 'Sediment module is run in offline mode' 211 IF (lwp) write(numsed,*) 'key_sed_off is activated at compilation time' 212 IF (lwp) write(numsed,*) 'ln_sed_2way is forced to false' 213 IF (lwp) write(numsed,*) '--------------------------------------------' 214 ln_sed_2way = .FALSE. 215 #endif 216 ! Initialisation of tracer damping 217 ! -------------------------------- 218 IF (ln_sediment_offline) THEN 219 CALL trc_dta_ini(jptra) 220 CALL trc_dmp_sed_ini 221 ENDIF 210 222 211 223 END SUBROUTINE sed_init 212 213 224 214 225 SUBROUTINE sed_init_geom … … 227 238 !! * Modules used 228 239 !! * 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 240 INTEGER :: ji, jj, jk, jn 241 REAL(wp) :: za0, za1, zt, zw, zsum, zsur, zprof, zprofw 242 REAL(wp) :: ztmp1, ztmp2 236 243 !---------------------------------------------------------- 237 244 238 WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry '239 WRITE(numsed,*) ' '245 IF(lwp) WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry ' 246 IF(lwp) WRITE(numsed,*) ' ' 240 247 241 248 ! Computation of 1D array of sediments points … … 250 257 END DO 251 258 259 IF ( indoce .EQ. 0 ) THEN 260 indoce = 1 261 iarroce(indoce) = 1 262 ENDIF 263 252 264 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 265 CALL ctl_stop( 'STOP', 'sed_ini: number of ocean points indoce doesn''t match number of point' ) 259 266 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 267 IF (lwp) WRITE(numsed,*) ' ' 268 IF (lwp) WRITE(numsed,*) ' total number of ocean points jpoce = ',jpoce 269 IF (lwp) WRITE(numsed,*) ' ' 270 ENDIF 285 271 286 272 ! initialization of dzkbot in [cm] … … 288 274 CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 289 275 dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 276 CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 290 277 291 278 ! Geometry and constants … … 293 280 ! (1st layer= diffusive layer = pur water) 294 281 !------------------------------------------ 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 282 za1 = ( sedzmin - sedhmax / FLOAT(jpksed-1) ) & 283 & / ( TANH((1-sedkth)/sedacr) - sedacr/FLOAT(jpksed-1) * ( LOG( COSH( (jpksed - sedkth) / sedacr) ) & 284 & - LOG( COSH( ( 1 - sedkth) / sedacr) ) ) ) 285 za0 = sedzmin - za1 * TANH( (1-sedkth) / sedacr ) 286 zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr ) ) 287 288 profsedw(1) = 0.0 289 profsed(1) = -dz(1) / 2. 290 DO jk = 2, jpksed 291 zw = REAL( jk , wp ) 292 zt = REAL( jk , wp ) - 0.5_wp 293 profsed(jk) = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) ) ) 294 profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) ) ) 295 END DO 296 297 dz(1) = 0.1 298 DO jk = 2, jpksed 299 dz(jk) = profsedw(jk) - profsedw(jk-1) 300 END DO 306 301 307 302 DO jk = 1, jpksed … … 311 306 ENDDO 312 307 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 308 ! Porosity profile [0] 321 309 !--------------------- 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 310 por(1) = 1.0 311 DO jk = 2, jpksed 312 por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk) ) ) 313 END DO 333 314 334 315 ! inverse of Porosity profile … … 353 334 dz3d(1:jpoce,1) = dz(1) 354 335 355 356 336 !--------------------------------------------- 357 337 ! Molecular weight [g/mol] for solid species 358 338 !--------------------------------------------- 359 360 339 361 340 ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16) … … 371 350 & 0.59 * 27. + 10. * 16. 372 351 352 mol_wgt(jsfeo) = 55.0 + 3.0 * ( 16.0 + 1.0) 353 354 mol_wgt(jsfes) = 55.0 + 32.0 355 373 356 ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1) 374 357 ! But den sity of Poc is an Hydrated material (= POC + 30H2O) … … 377 360 mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ & 378 361 & 16. * 14. + 31. ) / 122. 362 mol_wgt(jspos) = mol_wgt(jspoc) 363 mol_wgt(jspor) = mol_wgt(jspoc) 379 364 380 365 ! CaCO3 … … 384 369 ! Density of solid material in sediment [g/cm**3] 385 370 !------------------------------------------------ 386 dens = 2.6371 denssol = 2.6 387 372 388 373 ! Initialization of diffusion coefficient as function of porosity [cm**2/s] 389 374 !-------------------------------------------------------------------- 390 diff(:) = dcoef * por(:) 375 ! DO jn = 1, jpsol 376 ! DO jk = 1, jpksed 377 ! DO ji = 1, jpoce 378 ! diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) ) 379 ! END DO 380 ! END DO 381 ! END DO 382 383 ! Accumulation rate from Burwicz et al. (2011). This is used to 384 ! compute the flux of clays and minerals 385 ! -------------------------------------------------------------- 386 DO ji = 1, jpoce 387 ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 388 ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 ) 389 wacc(ji) = ztmp1 + ztmp2 390 END DO 391 391 392 392 393 393 ! Initialization of time step as function of porosity [cm**2/s] 394 394 !------------------------------------------------------------------ 395 rdtsed(2:jpksed) = dtsed / ( dens * por1(2:jpksed) )396 397 395 END SUBROUTINE sed_init_geom 398 396 … … 408 406 !!---------------------------------------------------------------------- 409 407 410 INTEGER :: & 411 numnamsed = 28 408 INTEGER :: numnamsed_ref = -1 !! Logical units for namelist sediment 409 INTEGER :: numnamsed_cfg = -1 !! Logical units for namelist sediment 410 INTEGER :: ios ! Local integer output status for namelist read 411 CHARACTER(LEN=20) :: clname 412 412 413 413 TYPE PSED … … 422 422 TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d 423 423 424 NAMELIST/nam_time/nfreq 424 NAMELIST/nam_run/nrseddt,ln_sed_2way 425 NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox 425 426 NAMELIST/nam_trased/sedsol, sedwat 426 427 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 428 NAMELIST/nam_inorg/rcopal, dcoef, rccal, ratligc, rcligc 429 NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs, & 430 & rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfes, rcfeso, & 431 & xksedo2, xksedno3, xksedfeo, xksedso4 432 NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, ln_irrig, xirrzsc 433 NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out 434 435 INTEGER :: ji, jn, jn1 436 436 !------------------------------------------------------- 437 437 438 WRITE(numsed,*) ' sed_init_nam : Read namelists '439 WRITE(numsed,*) ' '438 IF(lwp) WRITE(numsed,*) ' sed_init_nam : Read namelists ' 439 IF(lwp) WRITE(numsed,*) ' ' 440 440 441 441 ! ryear = 1 year converted in second 442 442 !------------------------------------ 443 WRITE(numsed,*) ' ' 444 WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear 445 WRITE(numsed,*) ' ' 443 IF (lwp) THEN 444 WRITE(numsed,*) ' ' 445 WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear 446 WRITE(numsed,*) ' ' 447 ENDIF 446 448 447 449 ! Reading namelist.sed variables 448 450 !--------------------------------- 449 CALL ctl_opn( numnamsed, 'namelist.sediment', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 450 451 dtsed = rdt 451 clname = 'namelist_sediment' 452 IF(lwp) WRITE(numsed,*) ' sed_init_nam : read SEDIMENT namelist' 453 IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~' 454 CALL ctl_opn( numnamsed_ref, TRIM( clname )//'_ref', 'OLD' , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 455 CALL ctl_opn( numnamsed_cfg, TRIM( clname )//'_cfg', 'OLD' , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 456 452 457 nitsed000 = nittrc000 453 458 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 ) 459 ! Namelist nam_run 460 REWIND( numnamsed_ref ) ! Namelist nam_run in reference namelist : Pisces variables 461 READ ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901) 462 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist', lwp ) 463 464 REWIND( numnamsed_cfg ) ! Namelist nam_run in reference namelist : Pisces variables 465 READ ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902) 466 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist', lwp ) 467 468 IF (lwp) THEN 469 WRITE(numsed,*) ' namelist nam_run' 470 WRITE(numsed,*) ' Nb of iterations for fast species nrseddt = ', nrseddt 471 WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way 472 ENDIF 473 474 REWIND( numnamsed_ref ) ! Namelist nam_geom in reference namelist : Pisces variables 475 READ ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903) 476 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist', lwp ) 477 478 REWIND( numnamsed_cfg ) ! Namelist nam_geom in reference namelist : Pisces variables 479 READ ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904) 480 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist', lwp ) 481 482 IF (lwp) THEN 483 WRITE(numsed,*) ' namelist nam_geom' 484 WRITE(numsed,*) ' Number of vertical layers jpksed = ', jpksed 485 WRITE(numsed,*) ' Minimum vertical spacing sedzmin = ', sedzmin 486 WRITE(numsed,*) ' Maximum depth of the sediment sedhmax = ', sedhmax 487 WRITE(numsed,*) ' Default parameter sedkth = ', sedkth 488 WRITE(numsed,*) ' Default parameter sedacr = ', sedacr 489 WRITE(numsed,*) ' Sediment porosity at the surface porsurf = ', porsurf 490 WRITE(numsed,*) ' Sediment porosity at infinite depth porinf = ', porinf 491 WRITE(numsed,*) ' Length scale of porosity variation rhox = ', rhox 492 ENDIF 493 494 jpksedm1 = jpksed - 1 495 dtsed = r2dttrc 496 497 REWIND( numnamsed_ref ) ! Namelist nam_trased in reference namelist : Pisces variables 498 READ ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905) 499 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist', lwp ) 500 501 REWIND( numnamsed_cfg ) ! Namelist nam_trased in reference namelist : Pisces variables 502 READ ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906) 503 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist', lwp ) 477 504 478 505 DO jn = 1, jpsol … … 489 516 END DO 490 517 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 518 IF (lwp) THEN 519 WRITE(numsed,*) ' namelist nam_trased' 520 WRITE(numsed,*) ' ' 521 DO jn = 1, jptrased 522 WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn)) 523 WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn)) 524 WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn)) 525 WRITE(numsed,*) ' ' 526 END DO 527 WRITE(numsed,*) ' ' 528 ENDIF 529 530 REWIND( numnamsed_ref ) ! Namelist nam_diased in reference namelist : Pisces variables 531 READ ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907) 532 907 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist', lwp ) 533 534 REWIND( numnamsed_cfg ) ! Namelist nam_diased in reference namelist : Pisces variables 535 READ ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908) 536 908 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist', lwp ) 501 537 502 REWIND( numnamsed )503 READ( numnamsed, nam_diased )504 505 538 DO jn = 1, jpdia3dsed 506 539 seddia3d(jn) = seddiag3d(jn)%snamesed … … 515 548 END DO 516 549 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 550 IF (lwp) THEN 551 WRITE(numsed,*) ' namelist nam_diased' 552 WRITE(numsed,*) ' ' 553 DO jn = 1, jpdia3dsed 554 WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn)) 555 WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn)) 556 WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn)) 557 WRITE(numsed,*) ' ' 558 END DO 559 560 DO jn = 1, jpdia2dsed 561 WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn)) 562 WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn)) 563 WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn)) 564 WRITE(numsed,*) ' ' 565 END DO 566 567 WRITE(numsed,*) ' ' 568 ENDIF 569 570 ! Inorganic chemistry parameters 537 571 !---------------------------------- 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 572 REWIND( numnamsed_ref ) ! Namelist nam_inorg in reference namelist : Pisces variables 573 READ ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909) 574 909 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist', lwp ) 575 576 REWIND( numnamsed_cfg ) ! Namelist nam_inorg in reference namelist : Pisces variables 577 READ ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910) 578 910 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist', lwp ) 579 580 IF (lwp) THEN 581 WRITE(numsed,*) ' namelist nam_inorg' 582 WRITE(numsed,*) ' reactivity for Si rcopal = ', rcopal 583 WRITE(numsed,*) ' diff. coef for por. dcoef = ', dcoef 584 WRITE(numsed,*) ' reactivity for calcite rccal = ', rccal 585 WRITE(numsed,*) ' L/C ratio in POC ratligc = ', ratligc 586 WRITE(numsed,*) ' reactivity for ligands rcligc = ', rcligc 587 WRITE(numsed,*) ' ' 588 ENDIF 548 589 549 590 ! Unity conversion to get saturation conc. psat in [mol.l-1] 550 591 ! and reactivity rc in [l.mol-1.s-1] 551 592 !---------------------------------------------------------- 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 593 reac_sil = rcopal / ryear 594 reac_ligc = rcligc / ryear 558 595 559 596 ! Additional parameter linked to POC/O2/No3/Po4 560 597 !---------------------------------------------- 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] 598 REWIND( numnamsed_ref ) ! Namelist nam_poc in reference namelist : Pisces variables 599 READ ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911) 600 911 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist', lwp ) 601 602 REWIND( numnamsed_cfg ) ! Namelist nam_poc in reference namelist : Pisces variables 603 READ ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912) 604 912 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist', lwp ) 605 606 IF (lwp) THEN 607 WRITE(numsed,*) ' namelist nam_poc' 608 WRITE(numsed,*) ' Redfield coef for oxy redO2 = ', redO2 609 WRITE(numsed,*) ' Redfield coef for no3 redNo3 = ', redNo3 610 WRITE(numsed,*) ' Redfield coef for po4 redPo4 = ', redPo4 611 WRITE(numsed,*) ' Redfield coef for carbon redC = ', redC 612 WRITE(numsed,*) ' Ration for iron bound P redfep = ', redfep 613 WRITE(numsed,*) ' reactivity for labile POC rcorgl = ', rcorgl 614 WRITE(numsed,*) ' reactivity for semi-refract. POC rcorgs = ', rcorgs 615 WRITE(numsed,*) ' reactivity for refractory POC rcorgr = ', rcorgr 616 WRITE(numsed,*) ' reactivity for NH4 rcnh4 = ', rcnh4 617 WRITE(numsed,*) ' reactivity for H2S rch2s = ', rch2s 618 WRITE(numsed,*) ' reactivity for Fe2+ rcfe2 = ', rcfe2 619 WRITE(numsed,*) ' reactivity for FeOH/H2S rcfeh2s = ', rcfeh2s 620 WRITE(numsed,*) ' reactivity for Fe2+/H2S rcfes = ', rcfes 621 WRITE(numsed,*) ' reactivity for FeS/O2 rcfeso = ', rcfeso 622 WRITE(numsed,*) ' Half-sat. cste for oxic remin xksedo2 = ', xksedo2 623 WRITE(numsed,*) ' Half-sat. cste for denit. xksedno3 = ', xksedno3 624 WRITE(numsed,*) ' Half-sat. cste for iron remin xksedfeo = ', xksedfeo 625 WRITE(numsed,*) ' Half-sat. cste for SO4 remin xksedso4 = ', xksedso4 626 WRITE(numsed,*) ' ' 627 ENDIF 628 629 630 so2ut = redO2 / redC 631 srno3 = redNo3 / redC 632 spo4r = redPo4 / redC 633 srDnit = ( (redO2 + 32. ) * 0.8 - redNo3 - redNo3 * 0.6 ) / redC 580 634 ! 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,*) ' ' 635 reac_pocl = rcorgl / ryear 636 reac_pocs = rcorgs / ryear 637 reac_pocr = rcorgr / ryear 638 reac_nh4 = rcnh4 / ryear 639 reac_h2s = rch2s / ryear 640 reac_fe2 = rcfe2 / ryear 641 reac_feh2s = rcfeh2s/ ryear 642 reac_fes = rcfes / ryear 643 reac_feso = rcfeso / ryear 591 644 592 645 ! reactivity rc in [l.mol-1.s-1] 593 646 reac_cal = rccal / ryear 594 647 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 648 ! Bioturbation parameter 607 649 !------------------------ 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. ) 650 REWIND( numnamsed_ref ) ! Namelist nam_btb in reference namelist : Pisces variables 651 READ ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913) 652 913 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist', lwp ) 653 654 REWIND( numnamsed_cfg ) ! Namelist nam_btb in reference namelist : Pisces variables 655 READ ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914) 656 914 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist', lwp ) 657 658 IF (lwp) THEN 659 WRITE(numsed,*) ' namelist nam_btb ' 660 WRITE(numsed,*) ' coefficient for bioturbation dbiot = ', dbiot 661 WRITE(numsed,*) ' Depth varying bioturbation ln_btbz = ', ln_btbz 662 WRITE(numsed,*) ' coefficient for btb attenuation dbtbzsc = ', dbtbzsc 663 WRITE(numsed,*) ' Adsorption coefficient of NH4 adsnh4 = ', adsnh4 664 WRITE(numsed,*) ' Bioirrigation in sediment ln_irrig = ', ln_irrig 665 WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc = ', xirrzsc 666 WRITE(numsed,*) ' ' 667 ENDIF 615 668 616 669 ! Initial value (t=0) for sediment pore water and solid components 617 670 !---------------------------------------------------------------- 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 ) 671 REWIND( numnamsed_ref ) ! Namelist nam_rst in reference namelist : Pisces variables 672 READ ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915) 673 915 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist', lwp ) 674 675 REWIND( numnamsed_cfg ) ! Namelist nam_rst in reference namelist : Pisces variables 676 READ ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916) 677 916 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist', lwp ) 678 679 IF (lwp) THEN 680 WRITE(numsed,*) ' namelist nam_rst ' 681 WRITE(numsed,*) ' boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 682 WRITE(numsed,*) ' ' 683 ENDIF 684 nn_dtsed = nn_dttrc 685 686 CLOSE( numnamsed_cfg ) 687 CLOSE( numnamsed_ref ) 624 688 625 689 END SUBROUTINE sed_init_nam 626 690 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 691 END MODULE sedini
Note: See TracChangeset
for help on using the changeset viewer.