Changeset 7938 for branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
- Timestamp:
- 2017-04-20T13:04:56+02:00 (8 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
r7927 r7938 98 98 USE bio_medusa_mod 99 99 USE bio_medusa_init_mod, ONLY: bio_medusa_init 100 USE carb_chem_mod, ONLY: carb_chem 101 USE bio_medusa_diag_slice_mod, ONLY: bio_medusa_diag_slice 100 102 USE bio_medusa_fin_mod, ONLY: bio_medusa_fin 101 USE bio_medusa_diag_slice_mod, ONLY: bio_medusa_diag_slice102 103 103 104 IMPLICIT NONE … … 160 161 !! 161 162 !! model state variables 162 REAL(wp), DIMENSION(jpi,jpj) :: zchn,zchd,zphn,zphd,zpds,zzmi 163 REAL(wp), DIMENSION(jpi,jpj) :: zzme,zdet,zdtc,zdin,zsil,zfer 164 REAL(wp) :: zage 165 # if defined key_roam 166 REAL(wp), DIMENSION(jpi,jpj) :: zdic, zalk, zoxy 167 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsal 168 # endif 169 # if defined key_mocsy 170 REAL(wp), DIMENSION(jpi,jpj) :: zpho 171 # endif 163 ! REAL(wp), DIMENSION(jpi,jpj) :: zchn,zchd,zphn,zphd,zpds,zzmi 164 ! REAL(wp), DIMENSION(jpi,jpj) :: zzme,zdet,zdtc,zdin,zsil,zfer 165 ! zage doesn't seem to be used - marc 19/4/17 166 ! REAL(wp) :: zage 167 !# if defined key_roam 168 ! REAL(wp), DIMENSION(jpi,jpj) :: zdic, zalk, zoxy 169 ! REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsal 170 !# endif 171 !# if defined key_mocsy 172 ! REAL(wp), DIMENSION(jpi,jpj) :: zpho 173 !# endif 172 174 !! 173 175 !! integrated source and sink terms … … 305 307 !! 306 308 !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have) 307 REAL(wp) :: f_xco2a308 REAL(wp), DIMENSION(jpi,jpj) :: f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux309 REAL(wp), DIMENSION(jpi,jpj) :: f_TDIC, f_TALK, f_dcf, f_henry310 REAL(wp), DIMENSION(jpi,jpj) :: f_pp0311 REAL(wp), DIMENSION(jpi,jpj) :: f_kw660, f_o2flux, f_o2sat309 ! REAL(wp) :: f_xco2a 310 ! REAL(wp), DIMENSION(jpi,jpj) :: f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux 311 ! REAL(wp), DIMENSION(jpi,jpj) :: f_TDIC, f_TALK, f_dcf, f_henry 312 ! REAL(wp), DIMENSION(jpi,jpj) :: f_pp0 313 ! REAL(wp), DIMENSION(jpi,jpj) :: f_kw660, f_o2flux, f_o2sat 312 314 REAL(wp) :: f_o2sat3 313 315 ! REAL(wp), DIMENSION(jpi,jpj) :: f_omcal, f_omarg 314 316 !! 315 317 !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen 316 REAL(wp), DIMENSION(jpi,jpj) :: f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm317 REAL(wp), DIMENSION(jpi,jpj) :: f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2318 ! REAL(wp), DIMENSION(jpi,jpj) :: f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm 319 ! REAL(wp), DIMENSION(jpi,jpj) :: f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2 318 320 !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s 319 321 REAL, PARAMETER :: weight_CO2_mol = 44.0095 !! g / mol … … 322 324 323 325 !! 324 INTEGER, DIMENSION(jpi,jpj) :: iters326 ! INTEGER, DIMENSION(jpi,jpj) :: iters 325 327 REAL(wp) :: f_year 326 328 INTEGER :: i_year … … 484 486 if (iyr1 .le. 1) then 485 487 !! before 1860 486 f_xco2a = hist_pco2(1)488 f_xco2a(:,:) = hist_pco2(1) 487 489 elseif (iyr2 .ge. 242) then 488 490 !! after 2099 489 f_xco2a = hist_pco2(242)491 f_xco2a(:,:) = hist_pco2(242) 490 492 else 491 493 !! just right … … 500 502 # endif 501 503 fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) 502 f_xco2a = fq4504 f_xco2a(:,:) = fq4 503 505 endif 504 506 # if defined key_axy_pi_co2 505 f_xco2a = 284.725 !! OCMIP pre-industrial pCO2507 f_xco2a(:,:) = 284.725 !! OCMIP pre-industrial pCO2 506 508 # endif 507 509 !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear =', nyear … … 513 515 !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2 =', fq2 514 516 !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3 =', fq3 515 IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2 =', f_xco2a 517 IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2 =', f_xco2a(1,1) 516 518 # endif 517 519 … … 544 546 !! field of omega calcite is used to determine the depth of the CCD 545 547 !!---------------------------------------------------------------------- 546 !! 547 IF(lwp) WRITE(numout,*) ' MEDUSA calculating all carbonate chemistry at kt =', kt 548 CALL flush(numout) 549 !! blank flags 550 i2_omcal(:,:) = 0 551 i2_omarg(:,:) = 0 552 !! loop over 3D space 553 DO jk = 1,jpk 554 DO jj = 2,jpjm1 555 DO ji = 2,jpim1 556 !! OPEN wet point IF..THEN loop 557 if (tmask(ji,jj,jk).eq.1) then 558 IF (lk_oasis) THEN 559 f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling 560 ENDIF 561 !! do carbonate chemistry 562 !! 563 !! set up required state variables 564 zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon 565 zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity 566 ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem) !! temperature 567 zsal(ji,jj) = tsn(ji,jj,jk,jp_sal) !! salinity 568 # if defined key_mocsy 569 zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) !! silicic acid 570 zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield 571 # endif 572 !! 573 !! AXY (28/02/14): check input fields 574 if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then 575 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 3D, ', & 576 tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', & 577 ji, ',', jj, ',', jk, ') at time', kt 578 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 3D, ', & 579 tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 580 ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) !! temperature 581 endif 582 if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then 583 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 3D, ', & 584 tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', & 585 ji, ',', jj, ',', jk, ') at time', kt 586 endif 587 !! 588 !! blank input variables not used at this stage (they relate to air-sea flux) 589 f_kw660(ji,jj) = 1.0 590 f_pp0(ji,jj) = 1.0 591 !! 592 !! calculate carbonate chemistry at grid cell midpoint 593 # if defined key_mocsy 594 !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 595 !! chemistry package 596 CALL mocsy_interface( ztmp(ji,jj), zsal(ji,jj), zalk(ji,jj), zdic(ji,jj), zsil(ji,jj), zpho(ji,jj), & ! inputs 597 f_pp0(ji,jj), fsdept(ji,jj,jk), gphit(ji,jj), f_kw660(ji,jj), f_xco2a, 1, & ! inputs 598 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 599 f_omcal(ji,jj), f_BetaD(ji,jj), f_rhosw(ji,jj), f_opres(ji,jj), f_insitut(ji,jj), & ! outputs 600 f_pco2atm(ji,jj), f_fco2atm(ji,jj), f_schmidtco2(ji,jj), f_kwco2(ji,jj), f_K0(ji,jj), & ! outputs 601 f_co2starair(ji,jj), f_co2flux(ji,jj), f_dpco2(ji,jj) ) ! outputs 602 !! 603 f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000. ! mmol / m3 -> umol / kg 604 f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000. ! meq / m3 -> ueq / kg 605 f_dcf(ji,jj) = f_rhosw(ji,jj) 606 # else 607 !! AXY (22/06/15): use old PML carbonate chemistry package (the 608 !! MEDUSA-2 default) 609 CALL trc_co2_medusa( ztmp(ji,jj), zsal(ji,jj), zdic(ji,jj), zalk(ji,jj), fsdept(ji,jj,jk), f_kw660(ji,jj), & ! inputs 610 f_xco2a, 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 611 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 612 !! 613 !! AXY (28/02/14): check output fields 614 if (iters(ji,jj) .eq. 25) then 615 IF(lwp) WRITE(numout,*) ' trc_bio_medusa: 3D ITERS WARNING, ', & 616 iters(ji,jj), ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 617 endif 618 # endif 619 !! 620 !! store 3D outputs 621 f3_pH(ji,jj,jk) = f_ph(ji,jj) 622 f3_h2co3(ji,jj,jk) = f_h2co3(ji,jj) 623 f3_hco3(ji,jj,jk) = f_hco3(ji,jj) 624 f3_co3(ji,jj,jk) = f_co3(ji,jj) 625 f3_omcal(ji,jj,jk) = f_omcal(ji,jj) 626 f3_omarg(ji,jj,jk) = f_omarg(ji,jj) 627 !! 628 !! CCD calculation: calcite 629 if (i2_omcal(ji,jj) .eq. 0 .and. f_omcal(ji,jj) .lt. 1.0) then 630 if (jk .eq. 1) then 631 f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk) 632 else 633 fq0 = f3_omcal(ji,jj,jk-1) - f_omcal(ji,jj) 634 fq1 = f3_omcal(ji,jj,jk-1) - 1.0 635 fq2 = fq1 / (fq0 + tiny(fq0)) 636 fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1) 637 fq4 = fq2 * fq3 638 f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4 639 endif 640 i2_omcal(ji,jj) = 1 641 endif 642 if ( i2_omcal(ji,jj) .eq. 0 .and. jk .eq. mbathy(ji,jj) ) then 643 !! reached seafloor and still no dissolution; set to seafloor (W-point) 644 f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1) 645 i2_omcal(ji,jj) = 1 646 endif 647 !! 648 !! CCD calculation: aragonite 649 if (i2_omarg(ji,jj) .eq. 0 .and. f_omarg(ji,jj) .lt. 1.0) then 650 if (jk .eq. 1) then 651 f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk) 652 else 653 fq0 = f3_omarg(ji,jj,jk-1) - f_omarg(ji,jj) 654 fq1 = f3_omarg(ji,jj,jk-1) - 1.0 655 fq2 = fq1 / (fq0 + tiny(fq0)) 656 fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1) 657 fq4 = fq2 * fq3 658 f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4 659 endif 660 i2_omarg(ji,jj) = 1 661 endif 662 if ( i2_omarg(ji,jj) .eq. 0 .and. jk .eq. mbathy(ji,jj) ) then 663 !! reached seafloor and still no dissolution; set to seafloor (W-point) 664 f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1) 665 i2_omarg(ji,jj) = 1 666 endif 667 endif 668 ENDDO 669 ENDDO 670 ENDDO 548 CALL carb_chem( kt ) 549 671 550 ENDIF 672 551 # endif … … 895 774 !!---------------------------------------------------------------------- 896 775 IF (lk_oasis) THEN 897 f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling776 f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling 898 777 ENDIF 899 778 !! … … 936 815 !! zero in this call 937 816 CALL mocsy_interface( ztmp(ji,jj), zsal(ji,jj), zalk(ji,jj), zdic(ji,jj), zsil(ji,jj), zpho(ji,jj), & ! inputs 938 f_pp0(ji,jj), 0.0, gphit(ji,jj), f_kw660(ji,jj), f_xco2a , 1, & ! inputs817 f_pp0(ji,jj), 0.0, gphit(ji,jj), f_kw660(ji,jj), f_xco2a(ji,jj), 1, & ! inputs 939 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 940 819 f_omcal(ji,jj), f_BetaD(ji,jj), f_rhosw(ji,jj), f_opres(ji,jj), f_insitut(ji,jj), & ! outputs … … 949 828 !! 950 829 !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not) 951 CALL trc_co2_medusa( ztmp(ji,jj), zsal(ji,jj), zdic(ji,jj), zalk(ji,jj), 0.0, f_kw660(ji,jj), f_xco2a , & ! inputs830 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 952 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 953 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 … … 968 847 !! quasi-sensible alternatives 969 848 f_ph(ji,jj) = 8.1 970 f_pco2w(ji,jj) = f_xco2a 849 f_pco2w(ji,jj) = f_xco2a(ji,jj) 971 850 f_h2co3(ji,jj) = 0.005 * zdic(ji,jj) 972 851 f_hco3(ji,jj) = 0.865 * zdic(ji,jj) … … 980 859 f_henry(ji,jj) = 1. 981 860 !! AXY (23/06/15): add in some extra MOCSY diagnostics 982 f_fco2w(ji,jj) = f_xco2a 861 f_fco2w(ji,jj) = f_xco2a(ji,jj) 983 862 f_BetaD(ji,jj) = 1. 984 863 f_rhosw(ji,jj) = 1.026 985 864 f_opres(ji,jj) = 0. 986 865 f_insitut(ji,jj) = ztmp(ji,jj) 987 f_pco2atm(ji,jj) = f_xco2a 988 f_fco2atm(ji,jj) = f_xco2a 866 f_pco2atm(ji,jj) = f_xco2a(ji,jj) 867 f_fco2atm(ji,jj) = f_xco2a(ji,jj) 989 868 f_schmidtco2(ji,jj) = 660. 990 869 f_kwco2(ji,jj) = 0. 991 870 f_K0(ji,jj) = 0. 992 f_co2starair(ji,jj) = f_xco2a 871 f_co2starair(ji,jj) = f_xco2a(ji,jj) 993 872 f_dpco2(ji,jj) = 0. 994 873 # endif … … 1110 989 !! AXY (24/11/16): add in extra MOCSY diagnostics 1111 990 IF( med_diag%ATM_XCO2%dgsave ) THEN 1112 f_xco2a_2d(ji,jj) = f_xco2a 991 f_xco2a_2d(ji,jj) = f_xco2a(ji,jj) 1113 992 ENDIF 1114 993 IF( med_diag%OCN_FCO2%dgsave ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
