MODULE air_sea_mod !!====================================================================== !! *** MODULE air_sea_mod *** !! Calculate the carbon chemistry for the whole ocean !!====================================================================== !! History : !! - ! 2017-04 (M. Stringer) Code taken from trcbio_medusa.F90 !!---------------------------------------------------------------------- #if defined key_medusa !!---------------------------------------------------------------------- !! MEDUSA bio-model !!---------------------------------------------------------------------- IMPLICIT NONE PRIVATE PUBLIC air_sea ! Called in trcbio_medusa.F90 !!---------------------------------------------------------------------- !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) !! $Id$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE air_sea( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE air_sea *** !! This called from TRC_BIO_MEDUSA and !! - calculate air-sea gas exchange !! - river inputs !!---------------------------------------------------------------------- USE bio_medusa_mod, ONLY: dms_andr, dms_andr2d, dms_aran, & dms_aran2d, dms_hall, dms_hall2d, & dms_simo, dms_simo2d, dms_surf, & dms_surf2d, iters, f_BetaD, & f_co2flux, f_co2flux2d, f_co2starair, & f_co2starair_2d, f_co3, & f_dcf, f_dpco2, f_fco2a_2d, f_fco2atm, & f_fco2w, f_fco2w_2d, f_h2co3, & f_hco3, f_henry, f_insitut, f_K0, & f_kw660, f_kw6602d, f_kwco2, f_kwo2, & f_o2flux, f_o2flux2d, & f_o2sat, f_o2sat2d, f_ocndpco2_2d, & f_ocnk0_2d, f_ocnkwco2_2d, & f_ocnrhosw_2d, f_ocnschco2_2d, & f_omarg, f_omcal, f_opres, & f_pco2a2d, f_pco2atm, & f_pco2w, f_pco2w2d, f_ph, f_pp0, & f_pp02d, f_rhosw, & f_riv_alk, f_riv_c, f_riv_n, & f_riv_si, f_runoff, & f_schmidtco2, f_TALK, f_TALK2d, & f_TDIC, f_TDIC2d, f_xco2a, f_xco2a_2d, & fgco2, & zalk, zchd, zchn, zdic, zdin, zoxy, & zpho, zsal, zsil, ztmp USE dom_oce, ONLY: e3t_0, e3t_n, gphit, tmask USE gastransfer, ONLY: gas_transfer USE iom, ONLY: lk_iomput USE in_out_manager, ONLY: lwp, numout USE mocsy_wrapper, ONLY: mocsy_interface USE oce, ONLY: PCO2a_in_cpl USE par_oce, ONLY: jpim1, jpjm1 USE sbc_oce, ONLY: fr_i, lk_oasis, qsr, wndm USE sms_medusa, ONLY: jdms, jdms_input, jdms_model, & jriver_alk, jriver_c, & jriver_n, jriver_si, & riv_alk, riv_c, riv_n, riv_si, & zn_dms_chd, zn_dms_chn, zn_dms_din, & zn_dms_mld, zn_dms_qsr USE trc, ONLY: med_diag USE trcco2_medusa, ONLY: trc_co2_medusa USE trcdms_medusa, ONLY: trc_dms_medusa USE trcoxy_medusa, ONLY: trc_oxy_medusa USE zdfmxl, ONLY: hmld !!* Substitution # include "domzgr_substitute.h90" !! time (integer timestep) INTEGER, INTENT( in ) :: kt !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s REAL, PARAMETER :: weight_CO2_mol = 44.0095 !! g / mol REAL, PARAMETER :: secs_in_day = 86400.0 !! s / d REAL, PARAMETER :: CO2flux_conv = (1.e-6 * weight_CO2_mol) / secs_in_day INTEGER :: ji, jj # if defined key_roam !!----------------------------------------------------------- !! Air-sea gas exchange !!----------------------------------------------------------- DO jj = 2,jpjm1 DO ji = 2,jpim1 !! OPEN wet point IF..THEN loop if (tmask(ji,jj,1) == 1) then IF (lk_oasis) THEN !! use 2D atm xCO2 from atm coupling f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj) ENDIF !! !! AXY (23/06/15): as part of an effort to update the !! carbonate chemistry in MEDUSA, the gas !! transfer velocity used in the carbon !! and oxygen cycles has been harmonised !! and is calculated by the same function !! here; this harmonisation includes !! changes to the PML carbonate chemistry !! scheme so that it too makes use of the !! same gas transfer velocity; the !! preferred parameterisation of this is !! Wanninkhof (2014), option 7 !! # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer' CALL flush(numout) # endif CALL gas_transfer( wndm(ji,jj), 1, 7, & ! inputs f_kw660(ji,jj) ) ! outputs # if defined key_debug_medusa IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' CALL flush(numout) # endif ENDIF ENDDO ENDDO DO jj = 2,jpjm1 DO ji = 2,jpim1 if (tmask(ji,jj,1) == 1) then !! air pressure (atm); ultimately this will use air !! pressure at the base of the UKESM1 atmosphere !! f_pp0(ji,jj) = 1.0 !! !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp =', ztmp(ji,jj) !! IF(lwp) WRITE(numout,*) ' MEDUSA wndm =', wndm(ji,jj) !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i =', fr_i(ji,jj) !! # if defined key_axy_carbchem # if defined key_mocsy !! !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate !! chemistry package; note that depth is set to !! zero in this call CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj), & zdic(ji,jj),zsil(ji,jj),zpho(ji,jj), & f_pp0(ji,jj),0.0, & gphit(ji,jj),f_kw660(ji,jj), & f_xco2a(ji,jj),1,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), & f_omcal(ji,jj),f_BetaD(ji,jj), & f_rhosw(ji,jj),f_opres(ji,jj), & f_insitut(ji,jj),f_pco2atm(ji,jj), & f_fco2atm(ji,jj),f_schmidtco2(ji,jj), & f_kwco2(ji,jj),f_K0(ji,jj), & f_co2starair(ji,jj),f_co2flux(ji,jj), & f_dpco2(ji,jj)) !! mmol / m3 -> umol / kg f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000. !! meq / m3 -> ueq / kg f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000. f_dcf(ji,jj) = f_rhosw(ji,jj) ENDIF ENDDO ENDDO # else DO jj = 2,jpjm1 DO ji = 2,jpim1 if (tmask(ji,jj,1) == 1) then iters(ji,jj) = 0 !! !! carbon dioxide (CO2); Jerry Blackford code (ostensibly !! OCMIP-2, but not) 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), & f_ph(ji,jj), & f_pco2w(ji,jj),f_h2co3(ji,jj), & f_hco3(ji,jj),f_co3(ji,jj), & f_omcal(ji,jj),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)) !! !! AXY (09/01/14): removed iteration and NaN checks; these have !! been moved to trc_co2_medusa together with a !! fudge that amends erroneous values (this is !! intended to be a temporary fudge!); the !! output warnings are retained here so that !! failure position can be determined if (iters(ji,jj) .eq. 25) then IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', & iters(ji,jj), ' AT (', ji, ', ', jj, ', 1) AT ', kt endif ENDIF ENDDO ENDDO # endif # else DO jj = 2,jpjm1 DO ji = 2,jpim1 if (tmask(ji,jj,1) == 1) then !! AXY (18/04/13): switch off carbonate chemistry !! calculations; provide quasi-sensible !! alternatives f_ph(ji,jj) = 8.1 f_pco2w(ji,jj) = f_xco2a(ji,jj) f_h2co3(ji,jj) = 0.005 * zdic(ji,jj) f_hco3(ji,jj) = 0.865 * zdic(ji,jj) f_co3(ji,jj) = 0.130 * zdic(ji,jj) f_omcal(ji,jj) = 4. f_omarg(ji,jj) = 2. f_co2flux(ji,jj) = 0. f_TDIC(ji,jj) = zdic(ji,jj) f_TALK(ji,jj) = zalk(ji,jj) f_dcf(ji,jj) = 1.026 f_henry(ji,jj) = 1. !! AXY (23/06/15): add in some extra MOCSY diagnostics f_fco2w(ji,jj) = f_xco2a(ji,jj) f_BetaD(ji,jj) = 1. f_rhosw(ji,jj) = 1.026 f_opres(ji,jj) = 0. f_insitut(ji,jj) = ztmp(ji,jj) f_pco2atm(ji,jj) = f_xco2a(ji,jj) f_fco2atm(ji,jj) = f_xco2a(ji,jj) f_schmidtco2(ji,jj) = 660. f_kwco2(ji,jj) = 0. f_K0(ji,jj) = 0. f_co2starair(ji,jj) = f_xco2a(ji,jj) f_dpco2(ji,jj) = 0. ENDIF ENDDO ENDDO # endif DO jj = 2,jpjm1 DO ji = 2,jpim1 if (tmask(ji,jj,1) == 1) then !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide !! through by layer thickness f_co2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_co2flux(ji,jj) * & 86400. / fse3t(ji,jj,1) !! !! oxygen (O2); OCMIP-2 code !! AXY (23/06/15): amend input list for oxygen to account !! for common gas transfer velocity !! Note that f_kwo2 is an about from the subroutine below, !! which doesn't seem to be used - marc 10/4/17 CALL trc_oxy_medusa(ztmp(ji,jj),zsal(ji,jj),f_kw660(ji,jj), & f_pp0(ji,jj),zoxy(ji,jj), & f_kwo2(ji,jj),f_o2flux(ji,jj), & f_o2sat(ji,jj)) !! !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide !! through by layer thickness f_o2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_o2flux(ji,jj) * & 86400. / fse3t(ji,jj,1) ENDIF ENDDO ENDDO !! Jpalm (08-2014) !! DMS surface concentration calculation !! initialy added for UKESM1 model. !! using MET-OFFICE subroutine. !! DMS module only needs Chl concentration and MLD !! to get an aproximate value of DMS concentration. !! air-sea fluxes are calculated by atmospheric chemitry model !! from atm and oc-surface concentrations. !! !! AXY (13/03/15): this is amended to calculate all of the DMS !! estimates examined during UKESM1 (see !! comments in trcdms_medusa.F90) !! IF (jdms == 1) THEN DO jj = 2,jpjm1 DO ji = 2,jpim1 if (tmask(ji,jj,1) == 1) then !! !! feed in correct inputs if (jdms_input .eq. 0) then !! use instantaneous inputs CALL trc_dms_medusa(zchn(ji,jj),zchd(ji,jj), & hmld(ji,jj),qsr(ji,jj), & zdin(ji,jj), & dms_andr(ji,jj),dms_simo(ji,jj), & dms_aran(ji,jj),dms_hall(ji,jj)) else !! use diel-average inputs CALL trc_dms_medusa(zn_dms_chn(ji,jj),zn_dms_chd(ji,jj), & zn_dms_mld(ji,jj),zn_dms_qsr(ji,jj), & zn_dms_din(ji,jj), & dms_andr(ji,jj),dms_simo(ji,jj), & dms_aran(ji,jj),dms_hall(ji,jj)) endif !! !! assign correct output to variable passed to atmosphere if (jdms_model .eq. 1) then dms_surf(ji,jj) = dms_andr(ji,jj) elseif (jdms_model .eq. 2) then dms_surf(ji,jj) = dms_simo(ji,jj) elseif (jdms_model .eq. 3) then dms_surf(ji,jj) = dms_aran(ji,jj) elseif (jdms_model .eq. 4) then dms_surf(ji,jj) = dms_hall(ji,jj) endif !! !! 2D diag through iom_use IF( lk_iomput ) THEN IF( med_diag%DMS_SURF%dgsave ) THEN dms_surf2d(ji,jj) = dms_surf(ji,jj) ENDIF IF( med_diag%DMS_ANDR%dgsave ) THEN dms_andr2d(ji,jj) = dms_andr(ji,jj) ENDIF IF( med_diag%DMS_SIMO%dgsave ) THEN dms_simo2d(ji,jj) = dms_simo(ji,jj) ENDIF IF( med_diag%DMS_ARAN%dgsave ) THEN dms_aran2d(ji,jj) = dms_aran(ji,jj) ENDIF IF( med_diag%DMS_HALL%dgsave ) THEN dms_hall2d(ji,jj) = dms_hall(ji,jj) ENDIF # if defined key_debug_medusa IF (lwp) write (numout,*) & 'trc_bio_medusa: finish calculating dms' CALL flush(numout) # endif ENDIF !! End iom ENDIF ENDDO ENDDO ENDIF !! End IF (jdms == 1) !! !! store 2D outputs !! !! JPALM -- 17-11-16 -- put fgco2 out of diag request !! is needed for coupling; pass through restart DO jj = 2,jpjm1 DO ji = 2,jpim1 if (tmask(ji,jj,1) == 1) then !! IF( med_diag%FGCO2%dgsave ) THEN !! convert from mol/m2/day to kg/m2/s !! mmol-C/m3/d -> kg-CO2/m2/s fgco2(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1) * & CO2flux_conv !! ENDIF IF ( lk_iomput ) THEN IF( med_diag%ATM_PCO2%dgsave ) THEN f_pco2a2d(ji,jj) = f_pco2atm(ji,jj) ENDIF IF( med_diag%OCN_PCO2%dgsave ) THEN f_pco2w2d(ji,jj) = f_pco2w(ji,jj) ENDIF IF( med_diag%CO2FLUX%dgsave ) THEN !! mmol/m3/d -> mmol/m2/d f_co2flux2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1) ENDIF IF( med_diag%TCO2%dgsave ) THEN f_TDIC2d(ji,jj) = f_TDIC(ji,jj) ENDIF IF( med_diag%TALK%dgsave ) THEN f_TALK2d(ji,jj) = f_TALK(ji,jj) ENDIF IF( med_diag%KW660%dgsave ) THEN f_kw6602d(ji,jj) = f_kw660(ji,jj) ENDIF IF( med_diag%ATM_PP0%dgsave ) THEN f_pp02d(ji,jj) = f_pp0(ji,jj) ENDIF IF( med_diag%O2FLUX%dgsave ) THEN f_o2flux2d(ji,jj) = f_o2flux(ji,jj) ENDIF IF( med_diag%O2SAT%dgsave ) THEN f_o2sat2d(ji,jj) = f_o2sat(ji,jj) ENDIF !! AXY (24/11/16): add in extra MOCSY diagnostics IF( med_diag%ATM_XCO2%dgsave ) THEN f_xco2a_2d(ji,jj) = f_xco2a(ji,jj) ENDIF IF( med_diag%OCN_FCO2%dgsave ) THEN f_fco2w_2d(ji,jj) = f_fco2w(ji,jj) ENDIF IF( med_diag%ATM_FCO2%dgsave ) THEN f_fco2a_2d(ji,jj) = f_fco2atm(ji,jj) ENDIF IF( med_diag%OCN_RHOSW%dgsave ) THEN f_ocnrhosw_2d(ji,jj) = f_rhosw(ji,jj) ENDIF IF( med_diag%OCN_SCHCO2%dgsave ) THEN f_ocnschco2_2d(ji,jj) = f_schmidtco2(ji,jj) ENDIF IF( med_diag%OCN_KWCO2%dgsave ) THEN f_ocnkwco2_2d(ji,jj) = f_kwco2(ji,jj) ENDIF IF( med_diag%OCN_K0%dgsave ) THEN f_ocnk0_2d(ji,jj) = f_K0(ji,jj) ENDIF IF( med_diag%CO2STARAIR%dgsave ) THEN f_co2starair_2d(ji,jj) = f_co2starair(ji,jj) ENDIF IF( med_diag%OCN_DPCO2%dgsave ) THEN f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj) ENDIF ENDIF ENDIF ENDDO ENDDO # endif !!----------------------------------------------------------------- !! River inputs !!----------------------------------------------------------------- DO jj = 2,jpjm1 DO ji = 2,jpim1 !! OPEN wet point IF..THEN loop if (tmask(ji,jj,1) == 1) then !! !! runoff comes in as kg / m2 / s !! used and written out as m3 / m2 / d (= m / d) !! where 1000 kg / m2 / d = !! 1 m3 / m2 / d = 1 m / d !! !! AXY (17/07/14): the compiler doesn't like this line for !! some reason; as MEDUSA doesn't even use !! runoff for riverine inputs, a temporary !! solution is to switch off runoff entirely !! here; again, this change is one of several !! that will need revisiting once MEDUSA has !! bedded down in UKESM1; particularly so if !! the land scheme provides information !! concerning nutrient fluxes !! !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * & !! 60. * 24. f_runoff(ji,jj) = 0.0 !! !! nutrients are added via rivers to the model in one of !! two ways: !! 1. via river concentration; i.e. the average nutrient !! concentration of a river water is described by a !! spatial file, and this is multiplied by runoff to !! give a nutrient flux !! 2. via direct river flux; i.e. the average nutrient !! flux due to rivers is described by a spatial file, !! and this is simply applied as a direct nutrient !! flux (i.e. it does not relate or respond to model !! runoff) nutrient fields are derived from the !! GlobalNEWS 2 database; carbon and alkalinity are !! derived from continent-scale DIC estimates (Huang et !! al., 2012) and some Arctic river alkalinity !! estimates (Katya?) !! !! as of 19/07/12, riverine nutrients can now be spread !! vertically across several grid cells rather than just !! poured into the surface box; this block of code is still !! executed, however, to set up the total amounts of !! nutrient entering via rivers !! !! nitrogen if (jriver_n .eq. 1) then !! river concentration specified; use runoff to !! calculate input f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj) elseif (jriver_n .eq. 2) then !! river flux specified; independent of runoff f_riv_n(ji,jj) = riv_n(ji,jj) endif !! !! silicon if (jriver_si .eq. 1) then !! river concentration specified; use runoff to !! calculate input f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj) elseif (jriver_si .eq. 2) then !! river flux specified; independent of runoff f_riv_si(ji,jj) = riv_si(ji,jj) endif !! !! carbon if (jriver_c .eq. 1) then !! river concentration specified; use runoff to !! calculate input f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj) elseif (jriver_c .eq. 2) then !! river flux specified; independent of runoff f_riv_c(ji,jj) = riv_c(ji,jj) endif !! !! alkalinity if (jriver_alk .eq. 1) then !! river concentration specified; use runoff to !! calculate input f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj) elseif (jriver_alk .eq. 2) then !! river flux specified; independent of runoff f_riv_alk(ji,jj) = riv_alk(ji,jj) endif ENDIF ENDDO ENDDO END SUBROUTINE air_sea #else !!====================================================================== !! Dummy module : No MEDUSA bio-model !!====================================================================== CONTAINS SUBROUTINE air_sea( ) ! Empty routine WRITE(*,*) 'air_sea: You should not have seen this print! error?' END SUBROUTINE air_sea #endif !!====================================================================== END MODULE air_sea_mod