Changeset 7958
- Timestamp:
- 2017-04-24T13:30:28+02:00 (7 years ago)
- Location:
- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA
- Files:
-
- 1 added
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
r7938 r7958 94 94 !! Oxygen production and consumption (and non-consumption) 95 95 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: foxy_prod,foxy_cons,foxy_anox 96 97 !! Add DMS in MEDUSA for UKESM1 model 98 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dms_surf 99 !! AXY (13/03/15): add in other DMS calculations 100 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dms_andr,dms_simo,dms_aran,dms_hall 96 101 97 102 !! Benthic fluxes … … 251 256 foxy_prod(jpi,jpj), foxy_cons(jpi,jpj), & 252 257 foxy_anox(jpi,jpj), & 258 dms_surf(jpi,jpj), & 259 dms_andr(jpi,jpj),dms_simo(jpi,jpj), & 260 dms_aran(jpi,jpj),dms_hall(jpi,jpj), & 253 261 f_sbenin_n(jpi,jpj),f_sbenin_fe(jpi,jpj), & 254 262 f_sbenin_c(jpi,jpj), & -
branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90
r7938 r7958 31 31 !! - ... 32 32 !!---------------------------------------------------------------------- 33 USE bio_medusa_mod, ONLY: iters, f_BetaD, f_dcf, f_fco2atm, & 34 f_co2flux, f_co2starair, f_co3, & 35 f_fco2w, f_dpco2, f_h2co3, f_hco3, & 33 USE bio_medusa_mod, ONLY: iters, f_BetaD, f_co2flux, & 34 f_co2starair, f_co3, f_dcf, f_dpco2, & 35 f_fco2atm, f_fco2w, f_h2co3, f_hco3, & 36 f_henry, & 36 37 f_insitut, f_K0, f_kw660, f_kwco2, & 37 38 f_omarg, f_omcal, f_opres, f_pco2atm, & … … 82 83 i2_omarg(:,:) = 0 83 84 84 !! loop over 3D space85 !! Loop over levels 85 86 DO jk = 1,jpk 87 86 88 DO jj = 2,jpjm1 87 89 DO ji = 2,jpim1 88 90 !! OPEN wet point IF..THEN loop 89 if (tmask(ji,jj,jk).eq.1) then91 IF (tmask(ji,jj,jk).eq.1) THEN 90 92 IF (lk_oasis) THEN 91 93 !! use 2D atm xCO2 from atm coupling … … 116 118 ' at (', ji, ',', jj, ',', jk, ') at time', kt 117 119 IF(lwp) WRITE(numout,*) & 118 ' carb_chem: T SWITCHING 3D, ', &120 ' carb_chem: T SWITCHING 3D, ', & 119 121 tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 120 122 ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) !! temperature … … 125 127 ' at (', ji, ',', jj, ',', jk, ') at time', kt 126 128 endif 127 ! MAYBE BUT A BREAK IN HERE - marc 128 !! 129 !! Blank input variables not used at this stage (they 130 !! relate to air-sea flux) 131 f_kw660(ji,jj) = 1.0 132 f_pp0(ji,jj) = 1.0 133 !! 129 ENDIF 130 ENDDO 131 ENDDO 132 133 !! Blank input variables not used at this stage (they 134 !! relate to air-sea flux) 135 f_kw660(:,:) = 1.0 136 f_pp0(:,:) = 1.0 137 138 DO jj = 2,jpjm1 139 DO ji = 2,jpim1 140 IF (tmask(ji,jj,jk).eq.1) THEN 134 141 !! calculate carbonate chemistry at grid cell midpoint 135 142 # if defined key_mocsy … … 168 175 f_omcal(ji,jj),f_omarg(ji,jj), & 169 176 f_co2flux(ji,jj),f_TDIC(ji,jj), & 170 f_TALK(ji,jj), f_dcf(ji,jj),&171 f_henry(ji,jj), 177 f_TALK(ji,jj),f_dcf(ji,jj), & 178 f_henry(ji,jj),iters(ji,jj)) 172 179 !! 173 180 !! AXY (28/02/14): check output fields … … 229 236 i2_omarg(ji,jj) = 1 230 237 endif 231 endif238 ENDIF 232 239 ENDDO 233 240 ENDDO -
branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
r7938 r7958 99 99 USE bio_medusa_init_mod, ONLY: bio_medusa_init 100 100 USE carb_chem_mod, ONLY: carb_chem 101 USE air_sea_mod, ONLY: air_sea 101 102 USE bio_medusa_diag_slice_mod, ONLY: bio_medusa_diag_slice 102 103 USE bio_medusa_fin_mod, ONLY: bio_medusa_fin … … 304 305 !! 305 306 !! flags to help with calculating the position of the CCD 306 INTEGER, DIMENSION(jpi,jpj) :: i2_omcal,i2_omarg 307 ! Moved into carb_chem.F90 - marc 20/4/17 308 ! INTEGER, DIMENSION(jpi,jpj) :: i2_omcal,i2_omarg 307 309 !! 308 310 !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have) … … 339 341 !! Jpalm (11-08-2014) 340 342 !! add DMS in MEDUSA for UKESM1 model 341 REAL(wp), DIMENSION(jpi,jpj) :: dms_surf343 ! REAL(wp), DIMENSION(jpi,jpj) :: dms_surf 342 344 !! AXY (13/03/15): add in other DMS calculations 343 REAL(wp), DIMENSION(jpi,jpj) :: dms_andr, dms_simo, dms_aran, dms_hall345 ! REAL(wp), DIMENSION(jpi,jpj) :: dms_andr, dms_simo, dms_aran, dms_hall 344 346 345 347 # endif … … 578 580 DO ji = 2,jpim1 579 581 !! OPEN wet point IF..THEN loop 580 if (tmask(ji,jj,jk) .eq.1) then581 !!=========================================================== ===========582 if (tmask(ji,jj,jk) == 1) then 583 !!=========================================================== 582 584 !! SETUP LOCAL GRID CELL 583 !!=========================================================== ===========584 !! 585 !!----------------------------------------------------------- ----------585 !!=========================================================== 586 !! 587 !!----------------------------------------------------------- 586 588 !! Some notes on grid vertical structure 587 !! - fsdepw(ji,jj,jk) is the depth of the upper surface of level jk 588 !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of level jk 589 !! - fsdepw(ji,jj,jk) is the depth of the upper surface of 590 !! level jk 591 !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of 592 !! level jk 589 593 !! - fse3t(ji,jj,jk) is the thickness of level jk 590 !!----------------------------------------------------------- ----------594 !!----------------------------------------------------------- 591 595 !! 592 596 !! AXY (01/03/10): set up level depth (bottom of level) 593 597 fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk) 594 598 !! AXY (28/11/16): local seafloor depth 595 !! previously mbathy(ji,jj) - 1, now mbathy(ji,jj) 599 !! previously mbathy(ji,jj) - 1, now 600 !! mbathy(ji,jj) 596 601 mbathy(ji,jj) = mbathy(ji,jj) 597 602 !! … … 599 604 !! negative values of state variables are not allowed to 600 605 !! contribute to the calculated fluxes 601 zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn)) !! non-diatom chlorophyll 602 zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd)) !! diatom chlorophyll 603 zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn)) !! non-diatoms 604 zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd)) !! diatoms 605 zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds)) !! diatom silicon 606 !! AXY (28/01/10): probably need to take account of chl/biomass connection 606 !! non-diatom chlorophyll 607 zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn)) 608 !! diatom chlorophyll 609 zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd)) 610 !! non-diatoms 611 zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn)) 612 !! diatoms 613 zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd)) 614 !! diatom silicon 615 zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds)) 616 !! AXY (28/01/10): probably need to take account of 617 !! chl/biomass connection 607 618 if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0. 608 619 if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0. … … 612 623 if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0. 613 624 if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0. 614 zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi)) !! microzooplankton 615 zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme)) !! mesozooplankton 616 zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet)) !! detrital nitrogen 617 zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) !! dissolved inorganic nitrogen 618 zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) !! dissolved silicic acid 619 zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer)) !! dissolved "iron" 625 !! microzooplankton 626 zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi)) 627 !! mesozooplankton 628 zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme)) 629 !! detrital nitrogen 630 zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet)) 631 !! dissolved inorganic nitrogen 632 zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) 633 !! dissolved silicic acid 634 zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) 635 !! dissolved "iron" 636 zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer)) 620 637 # if defined key_roam 621 zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) !! detrital carbon 622 zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon 623 zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity 624 zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) !! oxygen 638 !! detrital carbon 639 zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) 640 !! dissolved inorganic carbon 641 zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) 642 !! alkalinity 643 zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) 644 !! oxygen 645 zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) 625 646 # if defined key_axy_carbchem && defined key_mocsy 626 zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield 647 !! phosphate via DIN and Redfield 648 zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 627 649 # endif 628 650 !! … … 633 655 !! AXY (28/02/14): check input fields 634 656 if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then 635 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ', & 636 tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', & 637 ji, ',', jj, ',', jk, ') at time', kt 638 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ', & 639 tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 640 ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) !! temperature 657 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ', & 658 tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', & 659 ji, ',', jj, ',', jk, ') at time', kt 660 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ',& 661 tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 662 !! temperatur 663 ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) 641 664 endif 642 665 if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then 643 666 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', & 644 tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', &645 ji, ',', jj, ',', jk, ') at time', kt667 tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', & 668 ji, ',', jj, ',', jk, ') at time', kt 646 669 endif 647 670 # else 648 zdtc(ji,jj) = zdet(ji,jj) * xthetad !! implicit detrital carbon 671 !! implicit detrital carbon 672 zdtc(ji,jj) = zdet(ji,jj) * xthetad 649 673 # endif 650 674 # if defined key_debug_medusa 651 675 if (idf.eq.1) then 652 !! AXY (15/01/10)676 !! AXY (15/01/10) 653 677 if (trn(ji,jj,jk,jpdin).lt.0.) then 654 678 IF (lwp) write (numout,*) '------------------------------' 655 IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =', trn(ji,jj,jk,jpdin) 656 IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @', ji, jj, jk, kt 679 IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =', & 680 trn(ji,jj,jk,jpdin) 681 IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @', & 682 ji, jj, jk, kt 657 683 endif 658 684 if (trn(ji,jj,jk,jpsil).lt.0.) then 659 685 IF (lwp) write (numout,*) '------------------------------' 660 IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =', trn(ji,jj,jk,jpsil) 661 IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @', ji, jj, jk, kt 686 IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =', & 687 trn(ji,jj,jk,jpsil) 688 IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @', & 689 ji, jj, jk, kt 662 690 endif 663 691 # if defined key_roam 664 692 if (trn(ji,jj,jk,jpdic).lt.0.) then 665 693 IF (lwp) write (numout,*) '------------------------------' 666 IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =', trn(ji,jj,jk,jpdic) 667 IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @', ji, jj, jk, kt 694 IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =', & 695 trn(ji,jj,jk,jpdic) 696 IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @', & 697 ji, jj, jk, kt 668 698 endif 669 699 if (trn(ji,jj,jk,jpalk).lt.0.) then 670 700 IF (lwp) write (numout,*) '------------------------------' 671 IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =', trn(ji,jj,jk,jpalk) 672 IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @', ji, jj, jk, kt 701 IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =', & 702 trn(ji,jj,jk,jpalk) 703 IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @', & 704 ji, jj, jk, kt 673 705 endif 674 706 if (trn(ji,jj,jk,jpoxy).lt.0.) then 675 707 IF (lwp) write (numout,*) '------------------------------' 676 IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =', trn(ji,jj,jk,jpoxy) 677 IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @', ji, jj, jk, kt 708 IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =', & 709 trn(ji,jj,jk,jpoxy) 710 IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @', & 711 ji, jj, jk, kt 678 712 endif 679 713 # endif … … 713 747 IF( lk_iomput ) THEN 714 748 IF ( med_diag%INVTN%dgsave ) THEN 715 ftot_n(ji,jj) = ftot_n(ji,jj) + & 716 (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) + zzmi(ji,jj) + zzme(ji,jj) + zdet(ji,jj) + zdin(ji,jj) ) ) 749 ftot_n(ji,jj) = ftot_n(ji,jj) + & 750 (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) + & 751 zzmi(ji,jj) + zzme(ji,jj) + & 752 zdet(ji,jj) + zdin(ji,jj) ) ) 717 753 ENDIF 718 754 IF ( med_diag%INVTSI%dgsave ) THEN 719 ftot_si(ji,jj) = ftot_si(ji,jj) + &720 755 ftot_si(ji,jj) = ftot_si(ji,jj) + & 756 (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) ) 721 757 ENDIF 722 758 IF ( med_diag%INVTFE%dgsave ) THEN 723 ftot_fe(ji,jj) = ftot_fe(ji,jj) + & 724 (fse3t(ji,jj,jk) * ( xrfn * ( zphn(ji,jj) + zphd(ji,jj) + zzmi(ji,jj) + zzme(ji,jj) + zdet(ji,jj) ) + zfer(ji,jj) ) ) 759 ftot_fe(ji,jj) = ftot_fe(ji,jj) + & 760 (fse3t(ji,jj,jk) * ( xrfn * & 761 ( zphn(ji,jj) + zphd(ji,jj) + & 762 zzmi(ji,jj) + zzme(ji,jj) + & 763 zdet(ji,jj) ) + zfer(ji,jj) ) ) 725 764 ENDIF 726 765 # if defined key_roam 727 766 IF ( med_diag%INVTC%dgsave ) THEN 728 767 ftot_c(ji,jj) = ftot_c(ji,jj) + & 729 (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) + (xthetapd * zphd(ji,jj)) + & 730 (xthetazmi * zzmi(ji,jj)) + (xthetazme * zzme(ji,jj)) + zdtc(ji,jj) + & 731 zdic(ji,jj) ) ) 768 (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) + & 769 (xthetapd * zphd(ji,jj)) + & 770 (xthetazmi * zzmi(ji,jj)) + & 771 (xthetazme * zzme(ji,jj)) + & 772 zdtc(ji,jj) + zdic(ji,jj) ) ) 732 773 ENDIF 733 774 IF ( med_diag%INVTALK%dgsave ) THEN 734 ftot_a(ji,jj) = ftot_a(ji,jj) + (fse3t(ji,jj,jk) * ( zalk(ji,jj) ) ) 775 ftot_a(ji,jj) = ftot_a(ji,jj) + (fse3t(ji,jj,jk) * & 776 zalk(ji,jj)) 735 777 ENDIF 736 778 IF ( med_diag%INVTO2%dgsave ) THEN 737 ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) * ( zoxy(ji,jj) ) ) 779 ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) * & 780 zoxy(ji,jj)) 738 781 ENDIF 739 782 !! 740 783 !! AXY (10/11/16): CMIP6 diagnostics 741 784 IF ( med_diag%INTDISSIC%dgsave ) THEN 742 intdissic(ji,jj) = intdissic(ji,jj) + (fse3t(ji,jj,jk) * zdic(ji,jj)) 785 intdissic(ji,jj) = intdissic(ji,jj) + & 786 (fse3t(ji,jj,jk) * zdic(ji,jj)) 743 787 ENDIF 744 788 IF ( med_diag%INTDISSIN%dgsave ) THEN 745 intdissin(ji,jj) = intdissin(ji,jj) + (fse3t(ji,jj,jk) * zdin(ji,jj)) 789 intdissin(ji,jj) = intdissin(ji,jj) + & 790 (fse3t(ji,jj,jk) * zdin(ji,jj)) 746 791 ENDIF 747 792 IF ( med_diag%INTDISSISI%dgsave ) THEN 748 intdissisi(ji,jj) = intdissisi(ji,jj) + (fse3t(ji,jj,jk) * zsil(ji,jj)) 793 intdissisi(ji,jj) = intdissisi(ji,jj) + & 794 (fse3t(ji,jj,jk) * zsil(ji,jj)) 749 795 ENDIF 750 796 IF ( med_diag%INTTALK%dgsave ) THEN 751 inttalk(ji,jj) = inttalk(ji,jj) + (fse3t(ji,jj,jk) * zalk(ji,jj)) 797 inttalk(ji,jj) = inttalk(ji,jj) + & 798 (fse3t(ji,jj,jk) * zalk(ji,jj)) 752 799 ENDIF 753 800 IF ( med_diag%O2min%dgsave ) THEN … … 755 802 o2min(ji,jj) = zoxy(ji,jj) 756 803 IF ( med_diag%ZO2min%dgsave ) THEN 757 zo2min(ji,jj) = (fsdepw(ji,jj,jk) + fdep1(ji,jj)) / 2. !! layer midpoint 804 !! layer midpoint 805 zo2min(ji,jj) = (fsdepw(ji,jj,jk) + & 806 fdep1(ji,jj)) / 2.0 758 807 ENDIF 759 808 endif … … 764 813 CALL flush(numout) 765 814 766 !!====================================================================== 767 !! LOCAL GRID CELL CALCULATIONS 768 !!====================================================================== 769 !! 815 ENDIF 816 ENDDO 817 ENDDO 818 819 !!---------------------------------------------------------------- 820 !! Calculate air-sea gas exchange and river inputs 821 !!---------------------------------------------------------------- 822 IF ( jk == 1 ) THEN 823 call air_sea( kt ) 824 END IF 825 770 826 # if defined key_roam 771 if ( jk .eq. 1 ) then 772 !!---------------------------------------------------------------------- 773 !! Air-sea gas exchange 774 !!---------------------------------------------------------------------- 775 IF (lk_oasis) THEN 776 f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling 777 ENDIF 778 !! 779 !! AXY (23/06/15): as part of an effort to update the carbonate chemistry 780 !! in MEDUSA, the gas transfer velocity used in the carbon 781 !! and oxygen cycles has been harmonised and is calculated 782 !! by the same function here; this harmonisation includes 783 !! changes to the PML carbonate chemistry scheme so that 784 !! it too makes use of the same gas transfer velocity; the 785 !! preferred parameterisation of this is Wanninkhof (2014), 786 !! option 7 787 !! 788 # if defined key_debug_medusa 789 IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer' 790 CALL flush(numout) 791 # endif 792 CALL gas_transfer( wndm(ji,jj), 1, 7, & ! inputs 793 f_kw660(ji,jj) ) ! outputs 794 # if defined key_debug_medusa 795 IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' 796 CALL flush(numout) 797 # endif 798 !! 799 !! air pressure (atm); ultimately this will use air pressure at the base 800 !! of the UKESM1 atmosphere 801 !! 802 f_pp0(ji,jj) = 1.0 803 !! 804 !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp =', ztmp(ji,jj) 805 !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj) 806 !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_j =', zwind_j(ji,jj) 807 !! IF(lwp) WRITE(numout,*) ' MEDUSA wndm =', wndm(ji,jj) 808 !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i =', fr_i(ji,jj) 809 !! 810 # if defined key_axy_carbchem 811 # if defined key_mocsy 812 !! 813 !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 814 !! chemistry package; note that depth is set to 815 !! zero in this call 816 CALL mocsy_interface( ztmp(ji,jj), zsal(ji,jj), zalk(ji,jj), zdic(ji,jj), zsil(ji,jj), zpho(ji,jj), & ! inputs 817 f_pp0(ji,jj), 0.0, gphit(ji,jj), f_kw660(ji,jj), f_xco2a(ji,jj), 1, & ! inputs 818 f_ph(ji,jj), f_pco2w(ji,jj), f_fco2w(ji,jj), f_h2co3(ji,jj), f_hco3(ji,jj), f_co3(ji,jj), f_omarg(ji,jj), & ! outputs 819 f_omcal(ji,jj), f_BetaD(ji,jj), f_rhosw(ji,jj), f_opres(ji,jj), f_insitut(ji,jj), & ! outputs 820 f_pco2atm(ji,jj), f_fco2atm(ji,jj), f_schmidtco2(ji,jj), f_kwco2(ji,jj), f_K0(ji,jj), & ! outputs 821 f_co2starair(ji,jj), f_co2flux(ji,jj), f_dpco2(ji,jj) ) ! outputs 822 !! 823 f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000. ! mmol / m3 -> umol / kg 824 f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000. ! meq / m3 -> ueq / kg 825 f_dcf(ji,jj) = f_rhosw(ji,jj) 826 # else 827 iters(ji,jj) = 0 828 !! 829 !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not) 830 CALL trc_co2_medusa( ztmp(ji,jj), zsal(ji,jj), zdic(ji,jj), zalk(ji,jj), 0.0, f_kw660(ji,jj), f_xco2a(ji,jj), & ! inputs 831 f_ph(ji,jj), f_pco2w(ji,jj), f_h2co3(ji,jj), f_hco3(ji,jj), f_co3(ji,jj), f_omcal(ji,jj), & ! outputs 832 f_omarg(ji,jj), f_co2flux(ji,jj), f_TDIC(ji,jj), f_TALK(ji,jj), f_dcf(ji,jj), f_henry(ji,jj), iters(ji,jj) ) ! outputs 833 !! 834 !! AXY (09/01/14): removed iteration and NaN checks; these have 835 !! been moved to trc_co2_medusa together with a 836 !! fudge that amends erroneous values (this is 837 !! intended to be a temporary fudge!); the 838 !! output warnings are retained here so that 839 !! failure position can be determined 840 if (iters(ji,jj) .eq. 25) then 841 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', & 842 iters(ji,jj), ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 843 endif 844 # endif 845 # else 846 !! AXY (18/04/13): switch off carbonate chemistry calculations; provide 847 !! quasi-sensible alternatives 848 f_ph(ji,jj) = 8.1 849 f_pco2w(ji,jj) = f_xco2a(ji,jj) 850 f_h2co3(ji,jj) = 0.005 * zdic(ji,jj) 851 f_hco3(ji,jj) = 0.865 * zdic(ji,jj) 852 f_co3(ji,jj) = 0.130 * zdic(ji,jj) 853 f_omcal(ji,jj) = 4. 854 f_omarg(ji,jj) = 2. 855 f_co2flux(ji,jj) = 0. 856 f_TDIC(ji,jj) = zdic(ji,jj) 857 f_TALK(ji,jj) = zalk(ji,jj) 858 f_dcf(ji,jj) = 1.026 859 f_henry(ji,jj) = 1. 860 !! AXY (23/06/15): add in some extra MOCSY diagnostics 861 f_fco2w(ji,jj) = f_xco2a(ji,jj) 862 f_BetaD(ji,jj) = 1. 863 f_rhosw(ji,jj) = 1.026 864 f_opres(ji,jj) = 0. 865 f_insitut(ji,jj) = ztmp(ji,jj) 866 f_pco2atm(ji,jj) = f_xco2a(ji,jj) 867 f_fco2atm(ji,jj) = f_xco2a(ji,jj) 868 f_schmidtco2(ji,jj) = 660. 869 f_kwco2(ji,jj) = 0. 870 f_K0(ji,jj) = 0. 871 f_co2starair(ji,jj) = f_xco2a(ji,jj) 872 f_dpco2(ji,jj) = 0. 873 # endif 874 !! 875 !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide through by layer thickness 876 f_co2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_co2flux(ji,jj) * 86400. / fse3t(ji,jj,jk) 877 !! 878 !! oxygen (O2); OCMIP-2 code 879 !! AXY (23/06/15): amend input list for oxygen to account for common gas 880 !! transfer velocity 881 !! Note that f_kwo2 is an about from the subroutine below, 882 !! which doesn't seem to be used - marc 10/4/17 883 CALL trc_oxy_medusa( ztmp(ji,jj), zsal(ji,jj), f_kw660(ji,jj), f_pp0(ji,jj), zoxy(ji,jj), & ! inputs 884 f_kwo2(ji,jj), f_o2flux(ji,jj), f_o2sat(ji,jj) ) ! outputs 885 !! 886 !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide through by layer thickness 887 f_o2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_o2flux(ji,jj) * 86400. / fse3t(ji,jj,jk) 888 !! 889 !! Jpalm (08-2014) 890 !! DMS surface concentration calculation 891 !! initialy added for UKESM1 model. 892 !! using MET-OFFICE subroutine. 893 !! DMS module only needs Chl concentration and MLD 894 !! to get an aproximate value of DMS concentration. 895 !! air-sea fluxes are calculated by atmospheric chemitry model 896 !! from atm and oc-surface concentrations. 897 !! 898 !! AXY (13/03/15): this is amended to calculate all of the DMS 899 !! estimates examined during UKESM1 (see comments 900 !! in trcdms_medusa.F90) 901 !! 902 IF (jdms .eq. 1) THEN 903 !! 904 !! feed in correct inputs 905 if (jdms_input .eq. 0) then 906 !! use instantaneous inputs 907 CALL trc_dms_medusa( zchn(ji,jj), zchd(ji,jj), hmld(ji,jj), qsr(ji,jj), zdin(ji,jj), & ! inputs 908 dms_andr(ji,jj), dms_simo(ji,jj), dms_aran(ji,jj), dms_hall(ji,jj) ) ! outputs 909 else 910 !! use diel-average inputs 911 CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), & ! inputs 912 zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), zn_dms_din(ji,jj), & ! inputs 913 dms_andr(ji,jj), dms_simo(ji,jj), dms_aran(ji,jj), dms_hall(ji,jj) ) ! outputs 914 endif 915 !! 916 !! assign correct output to variable passed to atmosphere 917 if (jdms_model .eq. 1) then 918 dms_surf(ji,jj) = dms_andr(ji,jj) 919 elseif (jdms_model .eq. 2) then 920 dms_surf(ji,jj) = dms_simo(ji,jj) 921 elseif (jdms_model .eq. 3) then 922 dms_surf(ji,jj) = dms_aran(ji,jj) 923 elseif (jdms_model .eq. 4) then 924 dms_surf(ji,jj) = dms_hall(ji,jj) 925 endif 926 !! 927 !! 2D diag through iom_use 928 IF( lk_iomput ) THEN 929 IF( med_diag%DMS_SURF%dgsave ) THEN 930 dms_surf2d(ji,jj) = dms_surf(ji,jj) 931 ENDIF 932 IF( med_diag%DMS_ANDR%dgsave ) THEN 933 dms_andr2d(ji,jj) = dms_andr(ji,jj) 934 ENDIF 935 IF( med_diag%DMS_SIMO%dgsave ) THEN 936 dms_simo2d(ji,jj) = dms_simo(ji,jj) 937 ENDIF 938 IF( med_diag%DMS_ARAN%dgsave ) THEN 939 dms_aran2d(ji,jj) = dms_aran(ji,jj) 940 ENDIF 941 IF( med_diag%DMS_HALL%dgsave ) THEN 942 dms_hall2d(ji,jj) = dms_hall(ji,jj) 943 ENDIF 944 # if defined key_debug_medusa 945 IF (lwp) write (numout,*) 'trc_bio_medusa: finish calculating dms' 946 CALL flush(numout) 947 # endif 948 ENDIF 949 !! End iom 950 ENDIF 951 !! End DMS Loop 952 !! 953 !! store 2D outputs 954 !! 955 !! JPALM -- 17-11-16 -- put fgco2 out of diag request 956 !! is needed for coupling; pass through restart 957 !! IF( med_diag%FGCO2%dgsave ) THEN 958 !! convert from mol/m2/day to kg/m2/s 959 fgco2(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,jk) * CO2flux_conv !! mmol-C/m3/d -> kg-CO2/m2/s 960 !! ENDIF 961 IF ( lk_iomput ) THEN 962 IF( med_diag%ATM_PCO2%dgsave ) THEN 963 f_pco2a2d(ji,jj) = f_pco2atm(ji,jj) 964 ENDIF 965 IF( med_diag%OCN_PCO2%dgsave ) THEN 966 f_pco2w2d(ji,jj) = f_pco2w(ji,jj) 967 ENDIF 968 IF( med_diag%CO2FLUX%dgsave ) THEN 969 f_co2flux2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,jk) !! mmol/m3/d -> mmol/m2/d 970 ENDIF 971 IF( med_diag%TCO2%dgsave ) THEN 972 f_TDIC2d(ji,jj) = f_TDIC(ji,jj) 973 ENDIF 974 IF( med_diag%TALK%dgsave ) THEN 975 f_TALK2d(ji,jj) = f_TALK(ji,jj) 976 ENDIF 977 IF( med_diag%KW660%dgsave ) THEN 978 f_kw6602d(ji,jj) = f_kw660(ji,jj) 979 ENDIF 980 IF( med_diag%ATM_PP0%dgsave ) THEN 981 f_pp02d(ji,jj) = f_pp0(ji,jj) 982 ENDIF 983 IF( med_diag%O2FLUX%dgsave ) THEN 984 f_o2flux2d(ji,jj) = f_o2flux(ji,jj) 985 ENDIF 986 IF( med_diag%O2SAT%dgsave ) THEN 987 f_o2sat2d(ji,jj) = f_o2sat(ji,jj) 988 ENDIF 989 !! AXY (24/11/16): add in extra MOCSY diagnostics 990 IF( med_diag%ATM_XCO2%dgsave ) THEN 991 f_xco2a_2d(ji,jj) = f_xco2a(ji,jj) 992 ENDIF 993 IF( med_diag%OCN_FCO2%dgsave ) THEN 994 f_fco2w_2d(ji,jj) = f_fco2w(ji,jj) 995 ENDIF 996 IF( med_diag%ATM_FCO2%dgsave ) THEN 997 f_fco2a_2d(ji,jj) = f_fco2atm(ji,jj) 998 ENDIF 999 IF( med_diag%OCN_RHOSW%dgsave ) THEN 1000 f_ocnrhosw_2d(ji,jj) = f_rhosw(ji,jj) 1001 ENDIF 1002 IF( med_diag%OCN_SCHCO2%dgsave ) THEN 1003 f_ocnschco2_2d(ji,jj) = f_schmidtco2(ji,jj) 1004 ENDIF 1005 IF( med_diag%OCN_KWCO2%dgsave ) THEN 1006 f_ocnkwco2_2d(ji,jj) = f_kwco2(ji,jj) 1007 ENDIF 1008 IF( med_diag%OCN_K0%dgsave ) THEN 1009 f_ocnk0_2d(ji,jj) = f_K0(ji,jj) 1010 ENDIF 1011 IF( med_diag%CO2STARAIR%dgsave ) THEN 1012 f_co2starair_2d(ji,jj) = f_co2starair(ji,jj) 1013 ENDIF 1014 IF( med_diag%OCN_DPCO2%dgsave ) THEN 1015 f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj) 1016 ENDIF 1017 ENDIF 1018 !! 1019 endif 1020 !! End jk = 1 loop within ROAM key 827 828 ! Moved from above - marc 21/4/17 829 ! I think this should be moved into diagnostics at bottom - it 830 ! doesn't like it is used anyway else - marc 21/4/17 831 DO jj = 2,jpjm1 832 DO ji = 2,jpim1 833 !! OPEN wet point IF..THEN loop 834 if (tmask(ji,jj,jk) == 1) then 1021 835 1022 836 !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic … … 1025 839 o2sat3(ji, jj, jk) = f_o2sat3 1026 840 ENDIF 1027 1028 # endif 1029 1030 if ( jk .eq. 1 ) then 1031 !!---------------------------------------------------------------------- 1032 !! River inputs 1033 !!---------------------------------------------------------------------- 1034 !! 1035 !! runoff comes in as kg / m2 / s 1036 !! used and written out as m3 / m2 / d (= m / d) 1037 !! where 1000 kg / m2 / d = 1 m3 / m2 / d = 1 m / d 1038 !! 1039 !! AXY (17/07/14): the compiler doesn't like this line for some reason; 1040 !! as MEDUSA doesn't even use runoff for riverine inputs, 1041 !! a temporary solution is to switch off runoff entirely 1042 !! here; again, this change is one of several that will 1043 !! need revisiting once MEDUSA has bedded down in UKESM1; 1044 !! particularly so if the land scheme provides information 1045 !! concerning nutrient fluxes 1046 !! 1047 !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * 60. * 24. 1048 f_runoff(ji,jj) = 0.0 1049 !! 1050 !! nutrients are added via rivers to the model in one of two ways: 1051 !! 1. via river concentration; i.e. the average nutrient concentration 1052 !! of a river water is described by a spatial file, and this is 1053 !! multiplied by runoff to give a nutrient flux 1054 !! 2. via direct river flux; i.e. the average nutrient flux due to 1055 !! rivers is described by a spatial file, and this is simply applied 1056 !! as a direct nutrient flux (i.e. it does not relate or respond to 1057 !! model runoff) 1058 !! nutrient fields are derived from the GlobalNEWS 2 database; carbon and 1059 !! alkalinity are derived from continent-scale DIC estimates (Huang et al., 1060 !! 2012) and some Arctic river alkalinity estimates (Katya?) 1061 !! 1062 !! as of 19/07/12, riverine nutrients can now be spread vertically across 1063 !! several grid cells rather than just poured into the surface box; this 1064 !! block of code is still executed, however, to set up the total amounts 1065 !! of nutrient entering via rivers 1066 !! 1067 !! nitrogen 1068 if (jriver_n .eq. 1) then 1069 !! river concentration specified; use runoff to calculate input 1070 f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj) 1071 elseif (jriver_n .eq. 2) then 1072 !! river flux specified; independent of runoff 1073 f_riv_n(ji,jj) = riv_n(ji,jj) 1074 endif 1075 !! 1076 !! silicon 1077 if (jriver_si .eq. 1) then 1078 !! river concentration specified; use runoff to calculate input 1079 f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj) 1080 elseif (jriver_si .eq. 2) then 1081 !! river flux specified; independent of runoff 1082 f_riv_si(ji,jj) = riv_si(ji,jj) 1083 endif 1084 !! 1085 !! carbon 1086 if (jriver_c .eq. 1) then 1087 !! river concentration specified; use runoff to calculate input 1088 f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj) 1089 elseif (jriver_c .eq. 2) then 1090 !! river flux specified; independent of runoff 1091 f_riv_c(ji,jj) = riv_c(ji,jj) 1092 endif 1093 !! 1094 !! alkalinity 1095 if (jriver_alk .eq. 1) then 1096 !! river concentration specified; use runoff to calculate input 1097 f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj) 1098 elseif (jriver_alk .eq. 2) then 1099 !! river flux specified; independent of runoff 1100 f_riv_alk(ji,jj) = riv_alk(ji,jj) 1101 endif 1102 1103 endif 841 ENDIF 842 ENDDO 843 ENDDO 844 # endif 845 846 DO jj = 2,jpjm1 847 DO ji = 2,jpim1 848 !! OPEN wet point IF..THEN loop 849 if (tmask(ji,jj,jk) == 1) then 1104 850 1105 851 !!---------------------------------------------------------------------- … … 1370 1116 endif 1371 1117 # endif 1118 1119 ! MAYBE BUT A BREAK IN HERE, ZOOPLANKTON GRAZING - marc 20/4/17 1120 ! (plankton growth is 281 lines) 1372 1121 1373 1122 !!---------------------------------------------------------------------- … … 1520 1269 fmegrow(ji,jj) + (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) + fmeexcr(ji,jj) + ((1.0 - xbetan) * finme(ji,jj)) ) 1521 1270 1271 ! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PLANKTON LOSSES - marc 20/4/17 1272 ! (zoo plankton grazing is 152 lines) 1273 1522 1274 !!---------------------------------------------------------------------- 1523 1275 !! Plankton metabolic losses … … 1569 1321 if (jmzme.eq.4) fdzme(ji,jj) = xmzme * zzme(ji,jj) * & !! sigmoid 1570 1322 ((zzme(ji,jj) * zzme(ji,jj)) / (xkzme + (zzme(ji,jj) * zzme(ji,jj)))) 1323 1324 ! MAYBE BUT A BREAK IN HERE, DETRITUS PROCESS- marc 20/4/17 1325 ! (misc plankton loses is 53 lines - maybe we can link with another section) 1571 1326 1572 1327 !!---------------------------------------------------------------------- … … 1643 1398 endif 1644 1399 1400 ! MAYBE BUT A BREAK IN HERE, IRON CHEMISTRY AND SCAVENGING - marc 20/4/17 1401 ! (detritus processes is 74 lines) 1402 1645 1403 !!---------------------------------------------------------------------- 1646 1404 !! Iron chemistry and fractionation … … 1978 1736 # endif 1979 1737 1738 ! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PROCESSES - marc 20/4/17 1739 ! (iron chemistry and scavenging is 340 lines) 1740 1980 1741 !!---------------------------------------------------------------------- 1981 1742 !! Miscellaneous … … 2053 1814 ! (xthetazmi * fdzmi2(ji,jj)) + (xthetazme * fdzme2(ji,jj))) * fse3t(ji,jj,jk)) ! linear mortality 2054 1815 # endif 1816 1817 ! MAYBE BUT A BREAK IN HERE, FAST-SINKINIG DETRITUS - marc 20/4/17 1818 ! (miscellaneous processes is 79 lines) 2055 1819 2056 1820 !!---------------------------------------------------------------------- … … 2628 2392 CALL flush(numout) 2629 2393 2394 ! MAYBE BUT A BREAK IN HERE, BUSINESS AND UPDATING - marc 20/4/17 2395 ! (fast sinking detritus is 576 lines) 2396 2630 2397 !!====================================================================== 2631 2398 !! LOCAL GRID CELL TRENDS … … 3069 2836 !! endif 3070 2837 !! endif 2838 2839 ! MAYBE BUT A BREAK IN HERE - marc 20/4/17 2840 ! (this would make the previous section about 470 lines and the one below 2841 ! about 700 lines) 3071 2842 3072 2843 # if defined key_trc_diabio
Note: See TracChangeset
for help on using the changeset viewer.