Changeset 8012 for branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
- Timestamp:
- 2017-05-10T09:47:35+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
r7998 r8012 243 243 ! REAL(wp), DIMENSION(jpi,jpj) :: fslown, fslowc 244 244 ! REAL(wp), DIMENSION(jpi,jpj) :: fslownflux, fslowcflux 245 REAL(wp), DIMENSION(jpi,jpj) :: fregen,fregensi245 ! REAL(wp), DIMENSION(jpi,jpj) :: fregen,fregensi 246 246 ! REAL(wp), DIMENSION(jpi,jpj) :: fregenfast,fregenfastsi 247 247 # if defined key_roam … … 397 397 ENDIF 398 398 399 !!------------------------------------------------------------------ ----399 !!------------------------------------------------------------------ 400 400 !! b0 is present for debugging purposes; using b0 = 0 sets the tendency 401 401 !! terms of all biological equations to 0. 402 !!------------------------------------------------------------------ ----402 !!------------------------------------------------------------------ 403 403 !! 404 404 !! AXY (03/09/14): probably not the smartest move ever, but it'll fit … … 410 410 b0 = 1. 411 411 # endif 412 !!------------------------------------------------------------------ ----412 !!------------------------------------------------------------------ 413 413 !! fast detritus ballast scheme (0 = no; 1 = yes) 414 414 !! alternative to ballast scheme is same scheme but with no ballast 415 415 !! protection (not dissimilar to Martin et al., 1987) 416 !!------------------------------------------------------------------ ----416 !!------------------------------------------------------------------ 417 417 !! 418 418 iball = 1 419 419 420 !!------------------------------------------------------------------ ----420 !!------------------------------------------------------------------ 421 421 !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output 422 422 !! these should *only* be used in 1D since they give comprehensive 423 423 !! output for ecological functions in the model; primarily used in 424 424 !! debugging 425 !!------------------------------------------------------------------ ----425 !!------------------------------------------------------------------ 426 426 !! 427 427 idf = 0 … … 434 434 endif 435 435 436 !!------------------------------------------------------------------ ----436 !!------------------------------------------------------------------ 437 437 !! Initialise arrays to zero and set up arrays for diagnostics 438 !!---------------------------------------------------------------------- 439 ! tmp - marc 440 write(numout,*) 'bbb13. before call to bio_medusa_init,kt=',kt 441 flush(numout) 442 ! 438 !!------------------------------------------------------------------ 443 439 CALL bio_medusa_init( kt ) 444 ! tmp - marc 445 write(numout,*) 'bbb14. after call to bio_medusa_init,kt=',kt 446 flush(numout) 447 ! 448 !! 440 449 441 # if defined key_axy_nancheck 450 442 DO jn = 1,jptra … … 452 444 !! fq1 = MAXVAL(trn(:,:,:,jn)) 453 445 fq2 = SUM(trn(:,:,:,jn)) 454 !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK', & 455 !! & kt, jn, fq0, fq1, fq2 456 !! AXY (30/01/14): much to our surprise, the next line doesn't work on HECTOR 457 !! and has been replaced here with a specialist routine 446 !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK', & 447 !! kt, jn, fq0, fq1, fq2 448 !! AXY (30/01/14): much to our surprise, the next line doesn't 449 !! work on HECTOR and has been replaced here with 450 !! a specialist routine 458 451 !! if (fq2 /= fq2 ) then 459 452 if ( ieee_is_nan( fq2 ) ) then 460 453 !! there's a NaN here 461 if (lwp) write(numout,*) 'NAN detected in field', jn, 'at time', kt, 'at position:' 454 if (lwp) write(numout,*) 'NAN detected in field', jn, & 455 'at time', kt, 'at position:' 462 456 DO jk = 1,jpk 463 457 DO jj = 1,jpj … … 466 460 !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then 467 461 if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then 468 if (lwp) write (numout,'(a,1pe12.2,4i6)') 'NAN-CHECK',&469 &tmask(ji,jj,jk), ji, jj, jk, jn462 if (lwp) write (numout,'(a,1pe12.2,4i6)') & 463 'NAN-CHECK', tmask(ji,jj,jk), ji, jj, jk, jn 470 464 endif 471 465 enddo … … 479 473 480 474 # if defined key_debug_medusa 481 IF (lwp) write (numout,*) 'trc_bio_medusa: variables initialised and checked' 475 IF (lwp) write (numout,*) & 476 'trc_bio_medusa: variables initialised and checked' 482 477 CALL flush(numout) 483 478 # endif 484 479 485 480 # if defined key_roam 486 !!------------------------------------------------------------------ ----481 !!------------------------------------------------------------------ 487 482 !! calculate atmospheric pCO2 488 !!------------------------------------------------------------------ ----483 !!------------------------------------------------------------------ 489 484 !! 490 485 !! what's atmospheric pCO2 doing? (data start in 1859) … … 504 499 !! AXY (14/06/12): tweaked to make more sense (and be correct) 505 500 # if defined key_bs_axy_yrlen 506 fq3 = (real(nday_year) - 1.0 + fq2) / 360.0 !! bugfix: for 360d year with HadGEM2-ES forcing 501 !! bugfix: for 360d year with HadGEM2-ES forcing 502 fq3 = (real(nday_year) - 1.0 + fq2) / 360.0 507 503 # else 508 fq3 = (real(nday_year) - 1.0 + fq2) / 365.0 !! original use of 365 days (not accounting for leap year or 360d year) 504 !! original use of 365 days (not accounting for leap year or 505 !! 360d year) 506 fq3 = (real(nday_year) - 1.0 + fq2) / 365.0 509 507 # endif 510 508 fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) … … 512 510 endif 513 511 # if defined key_axy_pi_co2 514 f_xco2a(:,:) = 284.725 !! OCMIP pre-industrial pCO2 512 !! OCMIP pre-industrial pCO2 513 f_xco2a(:,:) = 284.725 515 514 # endif 516 515 !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear =', nyear … … 540 539 !!============================= 541 540 !! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call : 542 !! we don't want to call on the first time-step of all run submission, 543 !! but only on the very first time-step, and then every month 544 !! So we call on nittrc000 if not restarted run, 545 !! else if one month after last call. 546 !! assume one month is 30d --> 3600*24*30 : 2592000s 547 !! try to call carb-chem at 1st month's tm-stp : x * 30d + 1*rdt(i.e: mod = rdt) 541 !! we don't want to call on the first time-step of all run 542 !! submission, but only on the very first time-step, and 543 !! then every month. So we call on nittrc000 if not 544 !! restarted run, else if one month after last call. 545 !! Assume one month is 30d --> 3600*24*30 : 2592000s 546 !! try to call carb-chem at 1st month's tm-stp : 547 !! x * 30d + 1*rdt(i.e: mod = rdt) 548 548 !! ++ need to pass carb-chem output var through restarts 549 If ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR. mod(kt*rdt,2592000.) == rdt ) THEN 550 !!---------------------------------------------------------------------- 549 If ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR. & 550 ( mod(kt*rdt,2592000.) == rdt ) ) THEN 551 !!--------------------------------------------------------------- 551 552 !! Calculate the carbonate chemistry for the whole ocean on the first 552 553 !! simulation timestep and every month subsequently; the resulting 3D 553 554 !! field of omega calcite is used to determine the depth of the CCD 554 !!--------------------------------------------------------------- -------555 !!--------------------------------------------------------------- 555 556 CALL carb_chem( kt ) 556 557 … … 563 564 # endif 564 565 565 !!------------------------------------------------------------------ ----566 !!------------------------------------------------------------------ 566 567 !! MEDUSA has unified equation through the water column 567 568 !! (Diff. from LOBSTER which has two sets: bio- and non-bio layers) 568 569 !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1 569 !!------------------------------------------------------------------ ----570 !!------------------------------------------------------------------ 570 571 !! 571 572 !! NOTE: the ordering of the loops below differs from that of some other … … 583 584 !! OPEN horizontal loops 584 585 DO jj = 2,jpjm1 585 DO ji = 2,jpim1 586 !! OPEN wet point IF..THEN loop 587 if (tmask(ji,jj,jk) == 1) then 588 !!=========================================================== 589 !! SETUP LOCAL GRID CELL 590 !!=========================================================== 591 !! 592 !!----------------------------------------------------------- 593 !! Some notes on grid vertical structure 594 !! - fsdepw(ji,jj,jk) is the depth of the upper surface of 595 !! level jk 596 !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of 597 !! level jk 598 !! - fse3t(ji,jj,jk) is the thickness of level jk 599 !!----------------------------------------------------------- 600 !! 601 !! AXY (01/03/10): set up level depth (bottom of level) 602 fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk) 603 !! AXY (28/11/16): local seafloor depth 604 !! previously mbathy(ji,jj) - 1, now 605 !! mbathy(ji,jj) 606 ! I should be able to remove this - marc 5/5/17 607 ! mbathy(ji,jj) = mbathy(ji,jj) 608 !! 609 !! set up model tracers 610 !! negative values of state variables are not allowed to 611 !! contribute to the calculated fluxes 612 !! non-diatom chlorophyll 613 zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn)) 614 !! diatom chlorophyll 615 zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd)) 616 !! non-diatoms 617 zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn)) 618 !! diatoms 619 zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd)) 620 !! diatom silicon 621 zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds)) 622 !! AXY (28/01/10): probably need to take account of 623 !! chl/biomass connection 624 if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0. 625 if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0. 626 if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0. 627 if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0. 628 !! AXY (23/01/14): duh - why did I forget diatom silicon? 629 if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0. 630 if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0. 631 !! microzooplankton 632 zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi)) 633 !! mesozooplankton 634 zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme)) 635 !! detrital nitrogen 636 zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet)) 637 !! dissolved inorganic nitrogen 638 zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) 639 !! dissolved silicic acid 640 zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) 641 !! dissolved "iron" 642 zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer)) 586 DO ji = 2,jpim1 587 !! OPEN wet point IF..THEN loop 588 if (tmask(ji,jj,jk) == 1) then 589 !!========================================================= 590 !! SETUP LOCAL GRID CELL 591 !!========================================================= 592 !! 593 !!--------------------------------------------------------- 594 !! Some notes on grid vertical structure 595 !! - fsdepw(ji,jj,jk) is the depth of the upper surface of 596 !! level jk 597 !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of 598 !! level jk 599 !! - fse3t(ji,jj,jk) is the thickness of level jk 600 !!--------------------------------------------------------- 601 !! 602 !! AXY (01/03/10): set up level depth (bottom of level) 603 fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk) 604 !! 605 !! set up model tracers 606 !! negative values of state variables are not allowed to 607 !! contribute to the calculated fluxes 608 !! non-diatom chlorophyll 609 zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn)) 610 !! diatom chlorophyll 611 zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd)) 612 !! non-diatoms 613 zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn)) 614 !! diatoms 615 zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd)) 616 !! diatom silicon 617 zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds)) 618 !! AXY (28/01/10): probably need to take account of 619 !! chl/biomass connection 620 if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0. 621 if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0. 622 if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0. 623 if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0. 624 !! AXY (23/01/14): duh - why did I forget diatom silicon? 625 if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0. 626 if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0. 627 ENDIF 628 ENDDO 629 ENDDO 630 631 DO jj = 2,jpjm1 632 DO ji = 2,jpim1 633 if (tmask(ji,jj,1) == 1) then 634 !! microzooplankton 635 zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi)) 636 !! mesozooplankton 637 zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme)) 638 !! detrital nitrogen 639 zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet)) 640 !! dissolved inorganic nitrogen 641 zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) 642 !! dissolved silicic acid 643 zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) 644 !! dissolved "iron" 645 zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer)) 646 ENDIF 647 ENDDO 648 ENDDO 649 643 650 # if defined key_roam 644 !! detrital carbon 645 zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) 646 !! dissolved inorganic carbon 647 zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) 648 !! alkalinity 649 zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) 650 !! oxygen 651 zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) 651 DO jj = 2,jpjm1 652 DO ji = 2,jpim1 653 if (tmask(ji,jj,1) == 1) then 654 !! detrital carbon 655 zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) 656 !! dissolved inorganic carbon 657 zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) 658 !! alkalinity 659 zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) 660 !! oxygen 661 zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) 652 662 # if defined key_axy_carbchem && defined key_mocsy 653 !! phosphate via DIN and Redfield654 zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0663 !! phosphate via DIN and Redfield 664 zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 655 665 # endif 656 !! 657 !! also need physical parameters for gas exchange calculations 658 ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem) 659 zsal(ji,jj) = tsn(ji,jj,jk,jp_sal) 660 !! 661 !! AXY (28/02/14): check input fields 662 if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then 663 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ', & 664 tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', & 665 ji, ',', jj, ',', jk, ') at time', kt 666 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ',& 667 tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 668 !! temperatur 669 ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) 670 endif 671 if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then 672 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', & 673 tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', & 674 ji, ',', jj, ',', jk, ') at time', kt 675 endif 666 !! 667 !! also need physical parameters for gas exchange 668 !! calculations 669 ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem) 670 zsal(ji,jj) = tsn(ji,jj,jk,jp_sal) 671 !! 672 !! AXY (28/02/14): check input fields 673 if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then 674 IF(lwp) WRITE(numout,*) & 675 ' trc_bio_medusa: T WARNING 2D, ', & 676 tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), & 677 ' at (', ji, ',', jj, ',', jk, ') at time', kt 678 IF(lwp) WRITE(numout,*) & 679 ' trc_bio_medusa: T SWITCHING 2D, ', & 680 tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 681 !! temperatur 682 ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) 683 endif 684 if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then 685 IF(lwp) WRITE(numout,*) & 686 ' trc_bio_medusa: S WARNING 2D, ', & 687 tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), & 688 ' at (', ji, ',', jj, ',', jk, ') at time', kt 689 endif 690 ENDIF 691 ENDDO 692 ENDDO 676 693 # else 677 !! implicit detrital carbon 678 zdtc(ji,jj) = zdet(ji,jj) * xthetad 694 DO jj = 2,jpjm1 695 DO ji = 2,jpim1 696 if (tmask(ji,jj,1) == 1) then 697 !! implicit detrital carbon 698 zdtc(ji,jj) = zdet(ji,jj) * xthetad 699 ENDIF 700 ENDDO 701 ENDDO 679 702 # endif 680 703 # if defined key_debug_medusa 681 if (idf.eq.1) then 682 !! AXY (15/01/10) 683 if (trn(ji,jj,jk,jpdin).lt.0.) then 684 IF (lwp) write (numout,*) '------------------------------' 685 IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =', & 686 trn(ji,jj,jk,jpdin) 687 IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @', & 688 ji, jj, jk, kt 704 DO jj = 2,jpjm1 705 DO ji = 2,jpim1 706 if (tmask(ji,jj,1) == 1) then 707 if (idf.eq.1) then 708 !! AXY (15/01/10) 709 if (trn(ji,jj,jk,jpdin).lt.0.) then 710 IF (lwp) write (numout,*) & 711 '------------------------------' 712 IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =', & 713 trn(ji,jj,jk,jpdin) 714 IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @', & 715 ji, jj, jk, kt 716 endif 717 if (trn(ji,jj,jk,jpsil).lt.0.) then 718 IF (lwp) write (numout,*) & 719 '------------------------------' 720 IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =', & 721 trn(ji,jj,jk,jpsil) 722 IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @', & 723 ji, jj, jk, kt 724 endif 725 # if defined key_roam 726 if (trn(ji,jj,jk,jpdic).lt.0.) then 727 IF (lwp) write (numout,*) & 728 '------------------------------' 729 IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =', & 730 trn(ji,jj,jk,jpdic) 731 IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @', & 732 ji, jj, jk, kt 733 endif 734 if (trn(ji,jj,jk,jpalk).lt.0.) then 735 IF (lwp) write (numout,*) & 736 '------------------------------' 737 IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =', & 738 trn(ji,jj,jk,jpalk) 739 IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @', & 740 ji, jj, jk, kt 741 endif 742 if (trn(ji,jj,jk,jpoxy).lt.0.) then 743 IF (lwp) write (numout,*) & 744 '------------------------------' 745 IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =', & 746 trn(ji,jj,jk,jpoxy) 747 IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @', & 748 ji, jj, jk, kt 749 endif 750 # endif 689 751 endif 690 if (trn(ji,jj,jk,jpsil).lt.0.) then 691 IF (lwp) write (numout,*) '------------------------------' 692 IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =', & 693 trn(ji,jj,jk,jpsil) 694 IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @', & 695 ji, jj, jk, kt 696 endif 752 ENDIF 753 ENDDO 754 ENDDO 755 # endif 756 # if defined key_debug_medusa 757 ! I'M NOT SURE THIS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17 758 ! if (idf.eq.1.AND.idfval.eq.1) then 759 ! DO jj = 2,jpjm1 760 ! DO ji = 2,jpim1 761 ! if (tmask(ji,jj,1) == 1) then 762 ! !! report state variable values 763 ! IF (lwp) write (numout,*) & 764 ! '------------------------------' 765 ! IF (lwp) write (numout,*) 'fthk(',jk,') = ', & 766 ! fse3t(ji,jj,jk) 767 ! IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj) 768 ! IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj) 769 ! IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj) 770 ! IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj) 771 ! IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj) 772 ! IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj) 773 ! IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj) 774 ! IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj) 775 ! IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj) 697 776 # if defined key_roam 698 if (trn(ji,jj,jk,jpdic).lt.0.) then 699 IF (lwp) write (numout,*) '------------------------------' 700 IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =', & 701 trn(ji,jj,jk,jpdic) 702 IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @', & 703 ji, jj, jk, kt 704 endif 705 if (trn(ji,jj,jk,jpalk).lt.0.) then 706 IF (lwp) write (numout,*) '------------------------------' 707 IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =', & 708 trn(ji,jj,jk,jpalk) 709 IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @', & 710 ji, jj, jk, kt 711 endif 712 if (trn(ji,jj,jk,jpoxy).lt.0.) then 713 IF (lwp) write (numout,*) '------------------------------' 714 IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =', & 715 trn(ji,jj,jk,jpoxy) 716 IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @', & 717 ji, jj, jk, kt 718 endif 777 ! IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj) 778 ! IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj) 779 ! IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj) 780 ! IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj) 719 781 # endif 720 endif 721 # endif 782 ! ENDIF 783 ! ENDDO 784 ! ENDDO 785 ! endif 786 # endif 787 722 788 # if defined key_debug_medusa 723 !! report state variable values 724 if (idf.eq.1.AND.idfval.eq.1) then 725 IF (lwp) write (numout,*) '------------------------------' 726 IF (lwp) write (numout,*) 'fthk(',jk,') = ', fse3t(ji,jj,jk) 727 IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj) 728 IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj) 729 IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj) 730 IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj) 731 IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj) 732 IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj) 733 IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj) 734 IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj) 735 IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj) 736 # if defined key_roam 737 IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj) 738 IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj) 739 IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj) 740 IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj) 741 # endif 742 endif 743 # endif 744 745 # if defined key_debug_medusa 746 if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then 747 IF (lwp) write (numout,*) '------------------------------' 748 IF (lwp) write (numout,*) 'dust = ', dust(ji,jj) 749 endif 750 # endif 751 752 !! sum tracers for inventory checks 753 IF( lk_iomput ) THEN 754 IF ( med_diag%INVTN%dgsave ) THEN 755 ftot_n(ji,jj) = ftot_n(ji,jj) + & 756 (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) + & 757 zzmi(ji,jj) + zzme(ji,jj) + & 758 zdet(ji,jj) + zdin(ji,jj) ) ) 759 ENDIF 760 IF ( med_diag%INVTSI%dgsave ) THEN 761 ftot_si(ji,jj) = ftot_si(ji,jj) + & 762 (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) ) 763 ENDIF 764 IF ( med_diag%INVTFE%dgsave ) THEN 765 ftot_fe(ji,jj) = ftot_fe(ji,jj) + & 766 (fse3t(ji,jj,jk) * ( xrfn * & 767 ( zphn(ji,jj) + zphd(ji,jj) + & 768 zzmi(ji,jj) + zzme(ji,jj) + & 769 zdet(ji,jj) ) + zfer(ji,jj) ) ) 770 ENDIF 771 # if defined key_roam 772 IF ( med_diag%INVTC%dgsave ) THEN 773 ftot_c(ji,jj) = ftot_c(ji,jj) + & 774 (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) + & 775 (xthetapd * zphd(ji,jj)) + & 776 (xthetazmi * zzmi(ji,jj)) + & 777 (xthetazme * zzme(ji,jj)) + & 778 zdtc(ji,jj) + zdic(ji,jj) ) ) 779 ENDIF 780 IF ( med_diag%INVTALK%dgsave ) THEN 781 ftot_a(ji,jj) = ftot_a(ji,jj) + (fse3t(ji,jj,jk) * & 782 zalk(ji,jj)) 783 ENDIF 784 IF ( med_diag%INVTO2%dgsave ) THEN 785 ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) * & 786 zoxy(ji,jj)) 787 ENDIF 788 !! 789 !! AXY (10/11/16): CMIP6 diagnostics 790 IF ( med_diag%INTDISSIC%dgsave ) THEN 791 intdissic(ji,jj) = intdissic(ji,jj) + & 792 (fse3t(ji,jj,jk) * zdic(ji,jj)) 793 ENDIF 794 IF ( med_diag%INTDISSIN%dgsave ) THEN 795 intdissin(ji,jj) = intdissin(ji,jj) + & 796 (fse3t(ji,jj,jk) * zdin(ji,jj)) 797 ENDIF 798 IF ( med_diag%INTDISSISI%dgsave ) THEN 799 intdissisi(ji,jj) = intdissisi(ji,jj) + & 800 (fse3t(ji,jj,jk) * zsil(ji,jj)) 801 ENDIF 802 IF ( med_diag%INTTALK%dgsave ) THEN 803 inttalk(ji,jj) = inttalk(ji,jj) + & 804 (fse3t(ji,jj,jk) * zalk(ji,jj)) 805 ENDIF 806 IF ( med_diag%O2min%dgsave ) THEN 807 if ( zoxy(ji,jj) < o2min(ji,jj) ) then 808 o2min(ji,jj) = zoxy(ji,jj) 809 IF ( med_diag%ZO2min%dgsave ) THEN 810 !! layer midpoint 811 zo2min(ji,jj) = (fsdepw(ji,jj,jk) + & 812 fdep1(ji,jj)) / 2.0 813 ENDIF 814 endif 815 ENDIF 816 # endif 817 ENDIF 818 819 CALL flush(numout) 820 821 ENDIF 822 ENDDO 823 ENDDO 824 825 !!---------------------------------------------------------------- 789 ! I'M NOT SURE THIS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17 790 ! if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then 791 ! DO jj = 2,jpjm1 792 ! DO ji = 2,jpim1 793 ! if (tmask(ji,jj,1) == 1) then 794 ! IF (lwp) write (numout,*) & 795 ! '------------------------------' 796 ! IF (lwp) write (numout,*) 'dust = ', dust(ji,jj) 797 ! ENDIF 798 ! ENDDO 799 ! ENDDO 800 ! endif 801 # endif 802 803 !!--------------------------------------------------------------- 826 804 !! Calculate air-sea gas exchange and river inputs 827 !!--------------------------------------------------------------- -805 !!--------------------------------------------------------------- 828 806 IF ( jk == 1 ) THEN 829 call air_sea( kt ) 830 END IF 831 832 # if defined key_roam 833 834 ! Moved from above - marc 21/4/17 835 ! I think this should be moved into diagnostics at bottom - it 836 ! doesn't like it is used anyway else - marc 21/4/17 837 DO jj = 2,jpjm1 838 DO ji = 2,jpim1 839 !! OPEN wet point IF..THEN loop 840 if (tmask(ji,jj,jk) == 1) then 841 842 !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic 843 IF ( med_diag%O2SAT3%dgsave ) THEN 844 call oxy_sato( ztmp(ji,jj), zsal(ji,jj), f_o2sat3 ) 845 o2sat3(ji, jj, jk) = f_o2sat3 846 ENDIF 847 ENDIF 848 ENDDO 849 ENDDO 850 # endif 851 852 !!------------------------------------------------------------------ 807 CALL air_sea( kt ) 808 ENDIF 809 810 !!--------------------------------------------------------------- 853 811 !! Phytoplankton growth, zooplankton grazing and miscellaneous 854 812 !! plankton losses. 855 !!--------------------------------------------------------------- ---813 !!--------------------------------------------------------------- 856 814 CALL plankton( jk ) 857 815 858 !!--------------------------------------------------------------- -816 !!--------------------------------------------------------------- 859 817 !! Iron chemistry and scavenging 860 !!--------------------------------------------------------------- -818 !!--------------------------------------------------------------- 861 819 CALL iron_chem_scav( jk ) 862 820 863 ! Miscellaneous processes - marc 864 865 DO jj = 2,jpjm1 866 DO ji = 2,jpim1 867 !! OPEN wet point IF..THEN loop 868 if (tmask(ji,jj,jk) == 1) then 869 870 !!---------------------------------------------------------------------- 871 !! Miscellaneous 872 !!---------------------------------------------------------------------- 873 !! 874 !! diatom frustule dissolution 875 fsdiss(ji,jj) = xsdiss * zpds(ji,jj) 876 877 # if defined key_debug_medusa 878 !! report miscellaneous calculations 879 if (idf.eq.1.AND.idfval.eq.1) then 880 IF (lwp) write (numout,*) '------------------------------' 881 IF (lwp) write (numout,*) 'fsdiss(',jk,') = ', fsdiss(ji,jj) 882 endif 883 # endif 884 885 !!---------------------------------------------------------------------- 886 !! Slow detritus creation 887 !!---------------------------------------------------------------------- 888 !! this variable integrates the creation of slow sinking detritus 889 !! to allow the split between fast and slow detritus to be 890 !! diagnosed 891 fslown(ji,jj) = fdpn(ji,jj) + fdzmi(ji,jj) + ((1.0 - xfdfrac1) * fdpd(ji,jj)) + & 892 ((1.0 - xfdfrac2) * fdzme(ji,jj)) + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj))) 893 !! 894 !! this variable records the slow detrital sinking flux at this 895 !! particular depth; it is used in the output of this flux at 896 !! standard depths in the diagnostic outputs; needs to be 897 !! adjusted from per second to per day because of parameter vsed 898 fslownflux(ji,jj) = zdet(ji,jj) * vsed * 86400. 899 # if defined key_roam 900 !! 901 !! and the same for detrital carbon 902 fslowc(ji,jj) = (xthetapn * fdpn(ji,jj)) + (xthetazmi * fdzmi(ji,jj)) + & 903 (xthetapd * (1.0 - xfdfrac1) * fdpd(ji,jj)) + & 904 (xthetazme * (1.0 - xfdfrac2) * fdzme(ji,jj)) + & 905 ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj))) 906 !! 907 !! this variable records the slow detrital sinking flux at this 908 !! particular depth; it is used in the output of this flux at 909 !! standard depths in the diagnostic outputs; needs to be 910 !! adjusted from per second to per day because of parameter vsed 911 fslowcflux(ji,jj) = zdtc(ji,jj) * vsed * 86400. 912 # endif 913 914 !!---------------------------------------------------------------------- 915 !! Nutrient regeneration 916 !! this variable integrates total nitrogen regeneration down the 917 !! watercolumn; its value is stored and output as a 2D diagnostic; 918 !! the corresponding dissolution flux of silicon (from sources 919 !! other than fast detritus) is also integrated; note that, 920 !! confusingly, the linear loss terms from plankton compartments 921 !! are labelled as fdX2 when one might have expected fdX or fdX1 922 !!---------------------------------------------------------------------- 923 !! 924 !! nitrogen 925 fregen(ji,jj) = (( (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) + & ! messy feeding 926 (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) + & ! messy feeding 927 fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) + & ! excretion + D remin. 928 fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) + fdzme2(ji,jj)) * fse3t(ji,jj,jk)) ! linear mortality 929 !! 930 !! silicon 931 fregensi(ji,jj) = (( fsdiss(ji,jj) + ((1.0 - xfdfrac1) * fdpds(ji,jj)) + & ! dissolution + non-lin. mortality 932 ((1.0 - xfdfrac3) * fgmepds(ji,jj)) + & ! egestion by zooplankton 933 fdpds2(ji,jj)) * fse3t(ji,jj,jk)) ! linear mortality 934 # if defined key_roam 935 !! 936 !! carbon 937 ! Doesn't look this is used - marc 10/4/17 938 ! fregenc(ji,jj) = (( (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj))) + & ! messy feeding 939 ! (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) + & ! messy feeding 940 ! (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj))) + & ! messy feeding 941 ! fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) + & ! respiration + D remin. 942 ! (xthetapn * fdpn2(ji,jj)) + (xthetapd * fdpd2(ji,jj)) + & ! linear mortality 943 ! (xthetazmi * fdzmi2(ji,jj)) + (xthetazme * fdzme2(ji,jj))) * fse3t(ji,jj,jk)) ! linear mortality 944 # endif 945 946 ENDIF 947 ENDDO 948 ENDDO 949 950 !!------------------------------------------------------------------ 951 !! Detritus process 952 !!------------------------------------------------------------------ 821 !!--------------------------------------------------------------- 822 !! Detritus processes 823 !!--------------------------------------------------------------- 953 824 CALL detritus( jk, iball ) 954 825 955 !!--------------------------------------------------------------- ---826 !!--------------------------------------------------------------- 956 827 !! Updating tracers 957 !!--------------------------------------------------------------- ---828 !!--------------------------------------------------------------- 958 829 CALL bio_medusa_update( kt, jk ) 959 830 960 !!--------------------------------------------------------------- ---831 !!--------------------------------------------------------------- 961 832 !! Diagnostics 962 !!--------------------------------------------------------------- ---833 !!--------------------------------------------------------------- 963 834 CALL bio_medusa_diag( kt, jk ) 964 835 836 !!------------------------------------------------------- 837 !! 2d specific k level diags 838 !!------------------------------------------------------- 965 839 IF( lk_iomput .AND. .NOT. ln_diatrc ) THEN 966 967 !!-------------------------------------------------------968 !! 2d specific k level diags969 !!-------------------------------------------------------970 840 CALL bio_medusa_diag_slice( jk ) 971 972 841 ENDIF 973 842 … … 975 844 ENDDO 976 845 977 !!------------------------------------------------------------------ ----846 !!------------------------------------------------------------------ 978 847 !! Final calculations for diagnostics 979 !!------------------------------------------------------------------ ----848 !!------------------------------------------------------------------ 980 849 CALL bio_medusa_fin( kt ) 981 850
Note: See TracChangeset
for help on using the changeset viewer.
