New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7938 for branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 – NEMO

Ignore:
Timestamp:
2017-04-20T13:04:56+02:00 (8 years ago)
Author:
marc
Message:

Carbon chemistry has been pulled out of trcbio_medusa.F90 into carb_chem.F90

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  
    9898      USE bio_medusa_mod 
    9999      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 
    100102      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin 
    101       USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice 
    102103 
    103104      IMPLICIT NONE 
     
    160161      !! 
    161162      !! 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 
    172174      !! 
    173175      !! integrated source and sink terms 
     
    305307      !! 
    306308      !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have) 
    307       REAL(wp)                     ::    f_xco2a 
    308       REAL(wp), DIMENSION(jpi,jpj) ::    f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux 
    309       REAL(wp), DIMENSION(jpi,jpj) ::    f_TDIC, f_TALK, f_dcf, f_henry 
    310       REAL(wp), DIMENSION(jpi,jpj) ::    f_pp0 
    311       REAL(wp), DIMENSION(jpi,jpj) ::    f_kw660, f_o2flux, f_o2sat 
     309!      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 
    312314      REAL(wp)                     ::    f_o2sat3 
    313315!      REAL(wp), DIMENSION(jpi,jpj) ::    f_omcal, f_omarg 
    314316      !! 
    315317      !! 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_fco2atm 
    317       REAL(wp), DIMENSION(jpi,jpj) ::    f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2 
     318!      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 
    318320      !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s 
    319321      REAL, PARAMETER :: weight_CO2_mol = 44.0095  !! g / mol 
     
    322324 
    323325      !! 
    324       INTEGER, DIMENSION(jpi,jpj)  ::    iters 
     326!      INTEGER, DIMENSION(jpi,jpj)  ::    iters 
    325327      REAL(wp) ::    f_year 
    326328      INTEGER  ::    i_year 
     
    484486      if (iyr1 .le. 1) then 
    485487         !! before 1860 
    486          f_xco2a = hist_pco2(1) 
     488         f_xco2a(:,:) = hist_pco2(1) 
    487489      elseif (iyr2 .ge. 242) then 
    488490         !! after 2099 
    489          f_xco2a = hist_pco2(242) 
     491         f_xco2a(:,:) = hist_pco2(242) 
    490492      else 
    491493         !! just right 
     
    500502#  endif 
    501503         fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) 
    502          f_xco2a = fq4 
     504         f_xco2a(:,:) = fq4 
    503505      endif 
    504506#  if defined key_axy_pi_co2 
    505       f_xco2a = 284.725          !! OCMIP pre-industrial pCO2 
     507      f_xco2a(:,:) = 284.725          !! OCMIP pre-industrial pCO2 
    506508#  endif 
    507509      !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear 
     
    513515      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2 
    514516      !! 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) 
    516518# endif 
    517519 
     
    544546         !! field of omega calcite is used to determine the depth of the CCD 
    545547         !!---------------------------------------------------------------------- 
    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 
    671550      ENDIF 
    672551# endif 
     
    895774                  !!---------------------------------------------------------------------- 
    896775                  IF (lk_oasis) THEN 
    897                      f_xco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling 
     776                     f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling 
    898777                  ENDIF 
    899778                  !! 
     
    936815                  !!                 zero in this call 
    937816                  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,                   &  ! inputs 
     817                  f_pp0(ji,jj), 0.0, gphit(ji,jj), f_kw660(ji,jj), f_xco2a(ji,jj), 1,                   &  ! inputs 
    939818                  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 
    940819                  f_omcal(ji,jj), f_BetaD(ji,jj), f_rhosw(ji,jj), f_opres(ji,jj), f_insitut(ji,jj),            &  ! outputs 
     
    949828                  !! 
    950829                  !! 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,  &  ! inputs 
     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 
    952831                  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 
    953832                  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 
     
    968847                  !!                 quasi-sensible alternatives 
    969848                  f_ph(ji,jj)           = 8.1 
    970                   f_pco2w(ji,jj)        = f_xco2a 
     849                  f_pco2w(ji,jj)        = f_xco2a(ji,jj) 
    971850                  f_h2co3(ji,jj)        = 0.005 * zdic(ji,jj) 
    972851                  f_hco3(ji,jj)         = 0.865 * zdic(ji,jj) 
     
    980859                  f_henry(ji,jj)        = 1. 
    981860                  !! 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) 
    983862                  f_BetaD(ji,jj)        = 1. 
    984863                  f_rhosw(ji,jj)        = 1.026 
    985864                  f_opres(ji,jj)        = 0. 
    986865                  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) 
    989868                  f_schmidtco2(ji,jj)   = 660. 
    990869                  f_kwco2(ji,jj)        = 0. 
    991870                  f_K0(ji,jj)           = 0. 
    992                   f_co2starair(ji,jj)   = f_xco2a 
     871                  f_co2starair(ji,jj)   = f_xco2a(ji,jj) 
    993872                  f_dpco2(ji,jj)        = 0. 
    994873#  endif 
     
    1110989                      !! AXY (24/11/16): add in extra MOCSY diagnostics 
    1111990                      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) 
    1113992                      ENDIF 
    1114993                      IF( med_diag%OCN_FCO2%dgsave ) THEN 
Note: See TracChangeset for help on using the changeset viewer.