- Timestamp:
- 2017-12-06T17:10:46+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5518_GO6_conserv_Check/NEMOGCM/NEMO/TOP_SRC/trcini.F90
r8442 r8926 28 28 USE trcini_idtra ! idealize tracer initialisation 29 29 USE trcini_medusa ! MEDUSA initialisation 30 USE par_medusa ! MEDUSA parameters (needed for elemental cycles) 30 31 USE trcdta ! initialisation from files 31 32 USE daymod ! calendar manager … … 35 36 USE sbc_oce 36 37 USE trcice ! tracers in sea ice 37 38 USE sms_medusa ! MEDUSA initialisation 38 39 IMPLICIT NONE 39 40 PRIVATE … … 62 63 !! or read data or analytical formulation 63 64 !!--------------------------------------------------------------------- 64 INTEGER :: jk, jn, jl ! dummy loop indices 65 INTEGER :: ji, jj, jk, jn, jl ! dummy loop indices 66 # if defined key_medusa && defined key_roam 67 !! AXY (23/11/2017) 68 REAL(wp) :: zsum3d, zsum2d 69 REAL(wp) :: zq1, zq2, loc_vol, loc_area 70 REAL(wp), DIMENSION(6) :: loc_cycletot3, loc_cycletot2 71 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztot3d 72 REAL(wp), DIMENSION(jpi,jpj) :: ztot2d, carea 73 # endif 65 74 CHARACTER (len=25) :: charout 66 75 !!--------------------------------------------------------------------- … … 98 107 ! ! total volume of the ocean 99 108 areatot = glob_sum( cvol(:,:,:) ) 109 carea(:,:) = e1e2t(:,:) * tmask(:,:,1) 100 110 101 111 IF( lk_pisces ) CALL trc_ini_pisces ! PISCES bio-model … … 192 202 ENDIF 193 203 204 # if defined key_medusa && defined key_roam 205 ! AXY (17/11/2017): calculate initial totals of elemental cycles 206 ! 207 ! This is done in a very hard-wired way here; in future, this could be 208 ! replaced with loops and using a 2D array; one dimension would cover 209 ! the tracers, the other would be for the elements; each tracer would 210 ! have a factor for each element to say how much of that element was 211 ! in that tracer; for example, PHN would be 1.0 for N, xrfn for Fe and 212 ! xthetapn for C, with the other elements 0.0; the array entry for PHN 213 ! would then be (1. 0. xrfn xthetapn 0. 0.) for (N, Si, Fe, C, A, O2); 214 ! saving this for the next iteration 215 ! 216 cycletot(:) = 0._wp 217 ! report elemental totals at initialisation as we go along 218 IF ( lwp ) WRITE(numout,*) 219 IF ( lwp ) WRITE(numout,*) ' Elemental cycle totals: ' 220 IF ( lwp ) CALL flush(numout) 221 ! nitrogen 222 ztot3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 223 trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin) 224 ztot2d(:,:) = zn_sed_n(:,:) 225 zsum3d = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 226 zsum2d = glob_sum( ztot2d(:,:) * carea(:,:) ) 227 cycletot(1) = zsum3d + zsum2d 228 IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, cycletot(1) 229 IF ( lwp ) CALL flush(numout) 230 ! silicon 231 ztot3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil) 232 ztot2d(:,:) = zn_sed_si(:,:) 233 zsum3d = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 234 zsum2d = glob_sum( ztot2d(:,:) * carea(:,:) ) 235 cycletot(2) = zsum3d + zsum2d 236 IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, cycletot(2) 237 IF ( lwp ) CALL flush(numout) 238 ! iron 239 ztot3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 240 trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer) 241 ztot2d(:,:) = zn_sed_fe(:,:) 242 zsum3d = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 243 zsum2d = glob_sum( ztot2d(:,:) * carea(:,:) ) 244 cycletot(3) = zsum3d + zsum2d 245 IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, cycletot(3) 246 IF ( lwp ) CALL flush(numout) 247 ! carbon (uses fixed C:N ratios on plankton tracers) 248 ztot3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn) + (trn(:,:,:,jpphd) * xthetapd) + & 249 (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) + & 250 trn(:,:,:,jpdtc) + trn(:,:,:,jpdic) 251 ztot2d(:,:) = zn_sed_c(:,:) + zn_sed_ca(:,:) 252 zsum3d = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 253 zsum2d = glob_sum( ztot2d(:,:) * carea(:,:) ) 254 cycletot(4) = zsum3d + zsum2d 255 IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, cycletot(4) 256 IF ( lwp ) CALL flush(numout) 257 ! alkalinity (note benthic correction) 258 ztot3d(:,:,:) = trn(:,:,:,jpalk) 259 ztot2d(:,:) = zn_sed_ca(:,:) * 2._wp 260 zsum3d = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 261 zsum2d = glob_sum( ztot2d(:,:) * carea(:,:) ) 262 cycletot(5) = zsum3d + zsum2d 263 IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, cycletot(5) 264 IF ( lwp ) CALL flush(numout) 265 ! oxygen (note no benthic) 266 ztot3d(:,:,:) = trn(:,:,:,jpoxy) 267 ztot2d(:,:) = 0._wp 268 zsum3d = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 269 zsum2d = glob_sum( ztot2d(:,:) * carea(:,:) ) 270 cycletot(6) = zsum3d + zsum2d 271 IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, cycletot(6) 272 IF ( lwp ) CALL flush(numout) 273 ! Check 274 zsum3d = glob_sum( cvol(:,:,:) ) 275 zsum2d = glob_sum( carea(:,:) ) 276 IF ( lwp ) WRITE(numout,*) 277 IF ( lwp ) WRITE(numout,*) ' check : cvol : ', zsum3d 278 IF ( lwp ) WRITE(numout,*) ' check : carea : ', zsum2d 279 IF ( lwp ) WRITE(numout,*) 280 IF ( lwp ) CALL flush(numout) 281 ! 282 ! This second check dodges around the halo points in the grid 283 ! to check that glob_sum is doing what it's supposed to be 284 ! doing; note that the loop ordering here is the "correct" way 285 ! (according to Dr. Google) 286 ! 287 IF ( lwp ) WRITE(numout,*) 288 IF ( lwp ) WRITE(numout,*) ' Elemental cycle totals (check): ' 289 ! ocean cells 290 loc_cycletot3(:) = 0._wp 291 loc_vol = 0._wp 292 DO jk = 1, jpk 293 DO jj = nldj,nlej 294 DO ji = nldi,nlei 295 loc_vol = loc_vol + cvol(ji,jj,jk) 296 ! nitrogen 297 zq1 = trn(ji,jj,jk,jpphn) + trn(ji,jj,jk,jpphd) + trn(ji,jj,jk,jpzmi) + & 298 trn(ji,jj,jk,jpzme) + trn(ji,jj,jk,jpdet) + trn(ji,jj,jk,jpdin) 299 loc_cycletot3(1) = loc_cycletot3(1) + ( zq1 * cvol(ji,jj,jk) ) 300 ! silicon 301 zq1 = trn(ji,jj,jk,jppds) + trn(ji,jj,jk,jpsil) 302 loc_cycletot3(2) = loc_cycletot3(2) + ( zq1 * cvol(ji,jj,jk) ) 303 ! iron 304 zq1 = ((trn(ji,jj,jk,jpphn) + trn(ji,jj,jk,jpphd) + trn(ji,jj,jk,jpzmi) + & 305 trn(ji,jj,jk,jpzme) + trn(ji,jj,jk,jpdet)) * xrfn) + trn(ji,jj,jk,jpfer) 306 loc_cycletot3(3) = loc_cycletot3(3) + ( zq1 * cvol(ji,jj,jk) ) 307 ! carbon 308 zq1 = (trn(ji,jj,jk,jpphn) * xthetapn) + (trn(ji,jj,jk,jpphd) * xthetapd) + & 309 (trn(ji,jj,jk,jpzmi) * xthetazmi) + (trn(ji,jj,jk,jpzme) * xthetazme) + & 310 trn(ji,jj,jk,jpdtc) + trn(ji,jj,jk,jpdic) 311 loc_cycletot3(4) = loc_cycletot3(4) + ( zq1 * cvol(ji,jj,jk) ) 312 ! alkalinity 313 zq1 = trn(ji,jj,jk,jpalk) 314 loc_cycletot3(5) = loc_cycletot3(5) + ( zq1 * cvol(ji,jj,jk) ) 315 ! oxygen 316 zq1 = trn(ji,jj,jk,jpoxy) 317 loc_cycletot3(6) = loc_cycletot3(6) + ( zq1 * cvol(ji,jj,jk) ) 318 ENDDO 319 ENDDO 320 ENDDO 321 ! 322 ! sediment cells 323 loc_cycletot2(:) = 0._wp 324 loc_area = 0._wp 325 jk = 1 326 DO jj = nldj,nlej 327 DO ji = nldi,nlei 328 loc_area = loc_area + (e1e2t(ji,jj) * tmask(ji,jj,jk)) 329 ! nitrogen 330 zq1 = zn_sed_n(ji,jj) 331 loc_cycletot2(1) = loc_cycletot2(1) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 332 ! silicon 333 zq1 = zn_sed_si(ji,jj) 334 loc_cycletot2(2) = loc_cycletot2(2) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 335 ! iron 336 zq1 = zn_sed_fe(ji,jj) 337 loc_cycletot2(3) = loc_cycletot2(3) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 338 ! carbon 339 zq1 = zn_sed_c(ji,jj) + zn_sed_ca(ji,jj) 340 loc_cycletot2(4) = loc_cycletot2(4) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 341 ! alkalinity 342 zq1 = zn_sed_ca(ji,jj) * 2._wp 343 loc_cycletot2(5) = loc_cycletot2(5) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 344 ! skip oxygen 345 ENDDO 346 ENDDO 347 ! 348 ! report back 349 zq1 = loc_cycletot3(1) 350 zq2 = loc_cycletot2(1) 351 IF( lk_mpp ) CALL mpp_sum( zq1 ) 352 IF( lk_mpp ) CALL mpp_sum( zq2 ) 353 cycletot2(1) = zq1 + zq2 354 IF ( lwp ) WRITE(numout,9010) 'nitrogen', zq1, zq2, cycletot2(1) 355 zq1 = loc_cycletot3(2) 356 zq2 = loc_cycletot2(2) 357 IF( lk_mpp ) CALL mpp_sum( zq1 ) 358 IF( lk_mpp ) CALL mpp_sum( zq2 ) 359 cycletot2(2) = zq1 + zq2 360 IF ( lwp ) WRITE(numout,9010) 'silicon', zq1, zq2, cycletot2(2) 361 zq1 = loc_cycletot3(3) 362 zq2 = loc_cycletot2(3) 363 IF( lk_mpp ) CALL mpp_sum( zq1 ) 364 IF( lk_mpp ) CALL mpp_sum( zq2 ) 365 cycletot2(3) = zq1 + zq2 366 IF ( lwp ) WRITE(numout,9010) 'iron', zq1, zq2, cycletot2(3) 367 zq1 = loc_cycletot3(4) 368 zq2 = loc_cycletot2(4) 369 IF( lk_mpp ) CALL mpp_sum( zq1 ) 370 IF( lk_mpp ) CALL mpp_sum( zq2 ) 371 cycletot2(4) = zq1 + zq2 372 IF ( lwp ) WRITE(numout,9010) 'carbon', zq1, zq2, cycletot2(4) 373 zq1 = loc_cycletot3(5) 374 zq2 = loc_cycletot2(5) 375 IF( lk_mpp ) CALL mpp_sum( zq1 ) 376 IF( lk_mpp ) CALL mpp_sum( zq2 ) 377 cycletot2(5) = zq1 + zq2 378 IF ( lwp ) WRITE(numout,9010) 'alkalinity', zq1, zq2, cycletot2(5) 379 zq1 = loc_cycletot3(6) 380 zq2 = loc_cycletot2(6) 381 IF( lk_mpp ) CALL mpp_sum( zq1 ) 382 IF( lk_mpp ) CALL mpp_sum( zq2 ) 383 cycletot2(6) = zq1 + zq2 384 IF ( lwp ) WRITE(numout,9010) 'oxygen', zq1, zq2, cycletot2(6) 385 zq1 = loc_vol 386 zq2 = loc_area 387 IF( lk_mpp ) CALL mpp_sum( zq1 ) 388 IF( lk_mpp ) CALL mpp_sum( zq2 ) 389 IF ( lwp ) WRITE(numout,*) 390 IF ( lwp ) WRITE(numout,*) ' check : cvol : ', zq1 391 IF ( lwp ) WRITE(numout,*) ' check : carea : ', zq2 392 IF ( lwp ) WRITE(numout,*) 393 IF ( lwp ) CALL flush(numout) 394 ! 395 ! report back 396 ! nitrogen 397 ztot3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 398 trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin) 399 ztot2d(:,:) = zn_sed_n(:,:) 400 zsum3d = glob_sum( ztot3d(nldi:nlei,nldj:nlej,:) * cvol(nldi:nlei,nldj:nlej,:) ) 401 zsum2d = glob_sum( ztot2d(nldi:nlei,nldj:nlej) * carea(nldi:nlei,nldj:nlej) ) 402 cycletot(1) = zsum3d + zsum2d 403 IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, cycletot(1) 404 IF ( lwp ) CALL flush(numout) 405 406 ! 407 408 # endif 409 194 410 IF(lwp) WRITE(numout,*) 195 411 IF(lwp) WRITE(numout,*) 'trc_init : passive tracer set up completed' … … 202 418 203 419 9000 FORMAT(' tracer nb : ',i2,' name :',a10,' initial content :',e18.10) 420 9010 FORMAT(' element:',a10, & 421 ' 3d sum:',e18.10,' 2d sum:',e18.10, & 422 ' total:',e18.10) 204 423 ! 205 424 IF( nn_timing == 1 ) CALL timing_stop('trc_init')
Note: See TracChangeset
for help on using the changeset viewer.
