- Timestamp:
- 2021-02-03T16:03:34+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zsed.F90
r13546 r14385 66 66 REAL(wp) :: zsiloss, zcaloss, zws3, zws4, zwsc, zdep 67 67 REAL(wp) :: zwstpoc, zwstpon, zwstpop 68 REAL(wp) :: ztrfer, ztrpo4 s, ztrdp, zwdust, zmudia, ztemp68 REAL(wp) :: ztrfer, ztrpo4, ztrdop, zmudia, ztemp 69 69 REAL(wp) :: xdiano3, xdianh4 70 REAL(wp) :: zsoufer, zlight 70 71 ! 71 72 CHARACTER (len=25) :: charout 72 REAL(wp), DIMENSION(jpi,jpj ) :: zdenit2d, zbureff , zwork73 REAL(wp), DIMENSION(jpi,jpj ) :: zdenit2d, zbureff 73 74 REAL(wp), DIMENSION(jpi,jpj ) :: zwsbio3, zwsbio4 74 75 REAL(wp), DIMENSION(jpi,jpj ) :: zsedcal, zsedsi, zsedc 75 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight76 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep77 76 !!--------------------------------------------------------------------- 78 77 ! … … 81 80 82 81 ! Allocate temporary workspace 83 ALLOCATE( ztrpo4(jpi,jpj,jpk) )84 IF( ln_p5z ) ALLOCATE( ztrdop(jpi,jpj,jpk) )85 86 82 zdenit2d(:,:) = 0.e0 87 83 zbureff (:,:) = 0.e0 88 zwork (:,:) = 0.e089 84 zsedsi (:,:) = 0.e0 90 85 zsedcal (:,:) = 0.e0 91 86 zsedc (:,:) = 0.e0 92 87 88 ! OA: Warning, the following part is necessary to avoid CFL problems 89 ! above the sediments. Vertical sinking speed is limited using the 90 ! typical CFL criterion 91 ! -------------------------------------------------------------------- 92 DO_2D( 1, 1, 1, 1 ) 93 ikt = mbkt(ji,jj) 94 zdep = e3t(ji,jj,ikt,Kmm) / xstep 95 zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) ) 96 zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) ) 97 END_2D 98 93 99 IF( .NOT.lk_sed ) THEN 94 ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments 95 ! -------------------------------------------------------------------- 96 DO_2D( 1, 1, 1, 1 ) 97 ikt = mbkt(ji,jj) 98 zdep = e3t(ji,jj,ikt,Kmm) / xstep 99 zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) ) 100 zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) ) 101 END_2D 102 103 ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used 104 ! Computation of the fraction of organic matter that is permanently buried from Dunne's model 100 101 ! Computation of the sediment denitrification proportion: The metamodel 102 ! from Middleburg (2006) is used 103 ! Computation of the fraction of organic matter that is permanently 104 ! buried from Dunne's model (2007) 105 105 ! ------------------------------------------------------- 106 106 DO_2D( 1, 1, 1, 1 ) … … 125 125 ENDIF 126 126 127 ! This loss is scaled at each bottom grid cell for equilibrating the total budget of silica in the ocean.128 ! Thus, the amount of silica lost in the sediments equal the supply at the surface (dust+rivers)129 ! ------------------------------------------------------ 127 ! Fraction of dSi that is dissolved in the sediments. This fraction is 128 ! set to a constant value in p4zsbc 129 ! -------------------------------------------------------------------- 130 130 IF( .NOT.lk_sed ) zrivsil = 1._wp - sedsilfrac 131 131 132 ! Loss of bSi and CaCO3 to the sediments 132 133 DO_2D( 1, 1, 1, 1 ) 133 134 ikt = mbkt(ji,jj) … … 142 143 ! 143 144 IF( .NOT.lk_sed ) THEN 145 ! Dissolution of CaCO3 and bSi in the sediments. This is 146 ! instantaneous since here sediments are not explicitly 147 ! modeled. The amount of CaCO3 that dissolves in the sediments 148 ! is computed using a metamodel constructed from Archer (1996) 149 ! A minimum set to sedcalfrac is preserved. This value is defined in p4zbc 150 ! --------------------------------------------------------------- 144 151 DO_2D( 1, 1, 1, 1 ) 145 152 ikt = mbkt(ji,jj) … … 160 167 ENDIF 161 168 ! 169 ! Loss of particulate organic carbon and Fe to the sediments 162 170 DO_2D( 1, 1, 1, 1 ) 163 171 ikt = mbkt(ji,jj) … … 171 179 END_2D 172 180 ! 181 ! Loss of particulate organic N and P to the sediments (p5z) 173 182 IF( ln_p5z ) THEN 174 183 DO_2D( 1, 1, 1, 1 ) … … 185 194 186 195 IF( .NOT.lk_sed ) THEN 196 ! Degradation of organic matter in the sediments. The metamodel of 197 ! Middleburg (2006) is used here to mimic the diagenetic reactions. 187 198 ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after 188 199 ! denitrification in the sediments. Not very clever, but simpliest option. 200 ! ------------------------------------------------------------------------ 189 201 DO_2D( 1, 1, 1, 1 ) 190 202 ikt = mbkt(ji,jj) … … 192 204 zws4 = zwsbio4(ji,jj) * zdep 193 205 zws3 = zwsbio3(ji,jj) * zdep 206 ! Fraction that is permanently buried in the sediments 194 207 zrivno3 = 1. - zbureff(ji,jj) 195 208 zwstpoc = tr(ji,jj,ikt,jpgoc,Kbb) * zws4 + tr(ji,jj,ikt,jppoc,Kbb) * zws3 209 ! Denitrification in the sediments 196 210 zpdenit = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 211 ! Fraction that is not denitrified 197 212 z1pdenit = zwstpoc * zrivno3 - zpdenit 213 ! Oxic remineralization of organic matter in the sediments 198 214 zolimit = MIN( ( tr(ji,jj,ikt,jpoxy,Kbb) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 215 ! The fraction that cannot be denitrified nor oxidized by O2 216 ! is released back to the water column as DOC 199 217 tr(ji,jj,ikt,jpdoc,Krhs) = tr(ji,jj,ikt,jpdoc,Krhs) + z1pdenit - zolimit 218 ! Update of the tracers concentrations 200 219 tr(ji,jj,ikt,jppo4,Krhs) = tr(ji,jj,ikt,jppo4,Krhs) + zpdenit + zolimit 201 220 tr(ji,jj,ikt,jpnh4,Krhs) = tr(ji,jj,ikt,jpnh4,Krhs) + zpdenit + zolimit … … 206 225 sdenit(ji,jj) = rdenit * zpdenit * e3t(ji,jj,ikt,Kmm) 207 226 zsedc(ji,jj) = (1. - zrivno3) * zwstpoc * e3t(ji,jj,ikt,Kmm) 227 ! PISCES-QUOTA (p5z) 208 228 IF( ln_p5z ) THEN 209 229 zwstpop = tr(ji,jj,ikt,jpgop,Kbb) * zws4 + tr(ji,jj,ikt,jppop,Kbb) * zws3 … … 215 235 ENDIF 216 236 217 218 ! Nitrogen fixation process 219 ! Small source iron from particulate inorganic iron 220 !----------------------------------- 221 DO jk = 1, jpkm1 222 zlight (:,:,jk) = ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) ) 223 zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) ) 224 ENDDO 237 ! Nitrogen fixation process : light limitation of diazotrophy 238 ! Small source of iron from particulate inorganic iron (photochemistry) 239 ! This is a purely adhoc param. 240 !---------------------------------------------------------------------- 241 242 ! Diazotrophy (nitrogen fixation) is modeled according to an empirical 243 ! formulation. This is described in Aumont et al. (2015). Limitation 244 ! by P and Fe is computed. Inhibition by high N concentrations is imposed. 245 ! Diazotrophy sensitivity to temperature is parameterized as in 246 ! Ye et al. (2012) 247 ! ------------------------------------------------------------------------ 225 248 IF( ln_p4z ) THEN 249 ! PISCES part 226 250 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 227 ! !Potential nitrogen fixation dependant on temperature and iron251 ! Potential nitrogen fixation dependant on temperature and iron 228 252 ztemp = ts(ji,jj,jk,jp_tem,Kmm) 229 zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 230 ! Potential nitrogen fixation dependant on temperature and iron 253 zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) / rno3 254 ! Nitrogen fixation is inhibited when enough NO3 and/or NH4 255 zlim = ( 1.- xnanonh4(ji,jj,jk) - xnanono3(ji,jj,jk) ) 256 zfact = zlim * rfact2 257 ! Nitrogen fixation limitation by PO4 and Fe 258 ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 259 ztrpo4 = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 260 zlight = ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) * ( 1. - fr_i(ji,jj) ) 261 nitrpot(ji,jj,jk) = zmudia * r1_rday * zfact * MIN( ztrfer, ztrpo4 ) * zlight 262 END_3D 263 ELSE 264 ! PISCES-QUOTA part 265 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 266 ! Potential nitrogen fixation dependant on temperature and iron 267 ztemp = ts(ji,jj,jk,jp_tem,Kmm) 268 zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) / rno3 269 ! Nitrogen fixation is inhibited when enough NO3 and/or NH4 231 270 xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) ) 232 271 xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4) … … 234 273 IF( zlim <= 0.1 ) zlim = 0.01 235 274 zfact = zlim * rfact2 275 ! Nitrogen fixation limitation by PO4/DOP and Fe 236 276 ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 237 ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 238 ztrdp = ztrpo4(ji,jj,jk) 239 nitrpot(ji,jj,jk) = zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 240 END_3D 241 ELSE ! p5z 242 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 243 ! ! Potential nitrogen fixation dependant on temperature and iron 244 ztemp = ts(ji,jj,jk,jp_tem,Kmm) 245 zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 246 ! Potential nitrogen fixation dependant on temperature and iron 247 xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) ) 248 xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4) 249 zlim = ( 1.- xdiano3 - xdianh4 ) 250 IF( zlim <= 0.1 ) zlim = 0.01 251 zfact = zlim * rfact2 252 ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 253 ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 254 ztrdop(ji,jj,jk) = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4(ji,jj,jk)) 255 ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) 256 nitrpot(ji,jj,jk) = zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 277 ztrpo4 = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 278 ztrdop = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4) 279 zlight = ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) * ( 1. - fr_i(ji,jj) ) 280 nitrpot(ji,jj,jk) = zmudia * r1_rday * zfact * MIN( ztrfer, ztrpo4 + ztrdop ) * zlight 257 281 END_3D 258 282 ENDIF … … 261 285 ! ---------------------------------------- 262 286 IF( ln_p4z ) THEN 287 ! PISCES part 263 288 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 264 289 zfact = nitrpot(ji,jj,jk) * nitrfix 290 ! 1/3 of the diazotrophs growth is supposed to be excreted 291 ! as NH4. 1/3 as DOC and the rest is routed to POC/GOC as 292 ! a result of mortality by predation. Completely adhoc param 265 293 tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0 266 294 tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0 … … 270 298 tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0 271 299 tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 300 ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max 301 zlight = ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) * ( 1. - fr_i(ji,jj) ) 302 zsoufer = zlight * 2E-11 / ( 2E-11 + biron(ji,jj,jk) ) 272 303 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0 273 304 tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 274 305 tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 275 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer (ji,jj,jk)* rfact2 / rday306 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer * rfact2 / rday 276 307 tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + concdnh4 / ( concdnh4 + tr(ji,jj,jk,jppo4,Kbb) ) & 277 308 & * 0.001 * tr(ji,jj,jk,jpdoc,Kbb) * xstep 278 309 END_3D 279 310 ELSE ! p5z 311 ! PISCES-QUOTA part 280 312 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 281 313 zfact = nitrpot(ji,jj,jk) * nitrfix 314 ! 1/3 of the diazotrophs growth is supposed to be excreted 315 ! as NH4. 1/3 as DOC and the rest is routed POC and GOC as 316 ! a result of mortality by predation. Completely adhoc param 282 317 tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0 283 318 tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0 319 ! N/P ratio of diazotrophs is supposed to be 46 320 ztrpo4 = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 321 ztrdop = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4) 322 284 323 tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) & 285 & * ztrpo4 (ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk)+ rtrn)324 & * ztrpo4 / (ztrpo4 + ztrdop + rtrn) 286 325 tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zfact * 1.0 / 3.0 287 326 tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0 288 327 tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + 16.0 / 46.0 * zfact / 3.0 & 289 & - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk) & 290 & / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 328 & - 16.0 / 46.0 * zfact * ztrdop / (ztrpo4 + ztrdop + rtrn) 291 329 tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0 292 330 tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + zfact * 1.0 / 3.0 * 2.0 /3.0 … … 296 334 tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0 297 335 tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 336 ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max 337 zlight = ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) * ( 1. - fr_i(ji,jj) ) 338 zsoufer = zlight * 2E-11 / ( 2E-11 + biron(ji,jj,jk) ) 298 339 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0 299 340 tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 300 341 tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 301 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer (ji,jj,jk)* rfact2 / rday342 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer * rfact2 / rday 302 343 END_3D 303 344 ! … … 319 360 ENDIF 320 361 ! 321 IF( ln_p5z ) DEALLOCATE( ztrpo4, ztrdop )322 !323 362 IF( ln_timing ) CALL timing_stop('p4z_sed') 324 363 !
Note: See TracChangeset
for help on using the changeset viewer.
