!### macros ##################################################### #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if #include "tm5.inc" !################################################################# module emission_c13_ocean !----------------------------------------------------------------------- ! purpose ! ------- ! perform air-sea c13 emissions ! ! This version built for use with Takahashi et al. (2002) ! climatological delta pCO2, using Wanninkhof (1992) ! short-term quadratic gas transfer velocity and the Kettle et ! al. (2004) barometric correction. ! ! EMISSIONS ! interface: call emission_apply_ocean ! method: subroutine is called from its general parent routine apply_emission ! ! Table of contents: ! module emission_c13_ocean ! subroutine declare_emission_c13_ocean ! end subroutine declare_emission_c_ocean ! subroutine calc_emission_c13_ocean ! end subroutine calc_emission_c13_ocean ! subroutine emission_apply_c13_ocean(region) ! end subroutine emission_apply_c13_ocean ! subroutine free_emission_c13_ocean ! end subroutine free_emission_c13_ocean ! end module emission_c13_ocean ! !----------------------------------------------------------------------- use GO, only: gol, goErr, goPr, goLabel use Dims, only: nregions, nlon360,nlat180 use Meteo, only: u10m_dat,v10m_dat,ci_dat use Meteo, only: Set use chem_param, only: nmembersloc use global_types, only: emis_data use emission_common, only: enkf_is_forward, enkf_is_inverse, EmisInputDir use emission_common, only: co2_ocean, c13_ocean, c13_ocean_m, c13_ocean_n, flux_means_are_priors, fluxmultiplier use emission_data, only: do_add_2d, flux_to_gridbox_per_month implicit none ! --- const ------------------------- character(len=*), parameter :: mname = 'emission_c13_ocean' private ! public routines: public :: declare_emission_c13_ocean, free_emission_c13_ocean, calc_emission_c13_ocean, emission_apply_c13_ocean ! c13_ocean_regxmem is the result of the calculate routine. ! It holds the coarsened flux per ensemble member and region. type(emis_data),dimension(:,:),allocatable,target :: c13_ocean_regxmem character(len=200) :: deltapco2_prefix contains subroutine declare_emission_c13_ocean use GO, only: TrcFile, Init, Done, ReadRc use global_data, only: rcfile use dims, only: im,jm,nregions, iglbsfc implicit none integer :: region, imr, jmr, n integer :: status type(TrcFile) :: rcF character(len=*), parameter :: rname = mname//'/declare_emission_c13_ocean' if(enkf_is_forward()) then allocate(c13_ocean_regxmem(nregions,nmembersloc)) do region=1,nregions imr=im(region) jmr=jm(region) do n=1,nmembersloc allocate(c13_ocean_regxmem(region,n)%surf(imr,jmr)) end do enddo end if call Set( u10m_dat(iglbsfc), status, used=.true. ) call Set( v10m_dat(iglbsfc), status, used=.true. ) call Set( ci_dat(iglbsfc), status, used=.true. ) call Init( rcF, rcfile, status ) call ReadRc( rcF, 'deltapco2.prefix', deltapco2_prefix ,status) call Done( rcF ,status) end subroutine declare_emission_c13_ocean subroutine calc_emission_c13_ocean use Dims, only: dx,dy,nlon360, nlat180, im,jm, sec_month, newsrun, & idate, itau, staggered, nread, newmonth, newday, iglbsfc use datetime, only: calc_sm, tstamp, tau2date use emission_data, only: chardate, increase_itau use emission_common, only: ReadFromFile use toolbox, only: escape_tm use toolbox, only: coarsen_emission_2d use go_string, only: goNum2str use ParTools, only: myid, root ! for debugging use chem_param, only: xmc,xmc13 implicit none character(len=*), parameter :: rname = mname//'/calc_emission_c13_ocean' integer :: i,j,ii,jj,status integer, parameter :: annual_field = 0, rank2 = 2, level1 = 1 integer, dimension(6) :: idate_ice real,dimension(nlon360,nlat180),save :: ww,k92,schmidtno,solubility,deltapco2 real,pointer,dimension(:,:) :: u10m,v10m,ci integer, parameter :: add_field=0 real,dimension(nlon360,nlat180) :: c13_ocean_mem ! flux for a given ensemble member integer :: mem real,parameter :: delta_atm=-8. real :: alpha, epsilon character(len=10) :: cdate if(mod(itau,nread) /= 0) return ! only every nread hours ! get necessary fields if(newmonth) then cdate=trim(goNum2str(idate(1),'(i4.4)'))//trim(goNum2str(idate(2),'(i2.2)')) write(gol,'(a)') 'Ocean data to be read'//trim(cdate) ; call goPr call ReadFromFile(trim(EmisInputDir)//'/'//trim(deltapco2_prefix)//'.'//trim(goNum2str(idate(1)))//'.'//trim(goNum2str(idate(2),'(i2.2)'))//'.hdf', & 'DELTAPCO2', deltapco2, status) call ReadFromFile(trim(EmisInputDir)//'/'//trim(deltapco2_prefix)//'.'//trim(goNum2str(idate(1)))//'.'//trim(goNum2str(idate(2),'(i2.2)'))//'.hdf', & 'SOLUBILITY', solubility, status) call ReadFromFile(trim(EmisInputDir)//'/'//trim(deltapco2_prefix)//'.'//trim(goNum2str(idate(1)))//'.'//trim(goNum2str(idate(2),'(i2.2)'))//'.hdf', & 'SCHMIDTNO', schmidtno, status) if (status /= 0) call escape_tm('stop in '//rname) endif u10m=> u10m_dat(iglbsfc)%data(:,:,1) v10m=> v10m_dat(iglbsfc)%data(:,:,1) ci=> ci_dat(iglbsfc)%data(:,:,1) k92=0.0 ww=u10m**2+v10m**2 ! note: sqrt has been taken away since we square it in the next formula! ! make diffusion in m/s according to Wanninkhof92 param on 1x1 degree where(schmidtno(:,:).ge.0) k92=(0.31*(ww)/sqrt(schmidtno/660.0))*1.e-2/3600. ! only over oceans ! make fluxes of co2 in mol/m2/s with tm5 meteo. 1e-6 to change deltapco2 from uatm to atm. co2_ocean=0 where((solubility(:,:).ge.0).and.(deltapco2(:,:).gt.-1e10)) co2_ocean=k92*solubility*deltapco2*1e-6 co2_ocean=co2_ocean*(1.0-ci) ! multiply out the sea-ice for this day epsilon=-2. alpha=1+epsilon/1000 c13_ocean = co2_ocean*(0.011112*(1.+delta_atm/1000.))*alpha ! accumulate weekly average if(allocated(c13_ocean_m)) then if(flux_means_are_priors .eqv. .TRUE.) then c13_ocean_m=c13_ocean_m+c13_ocean ! mol/m2/s else c13_ocean_m=c13_ocean_m+c13_ocean*fluxmultiplier(:,:,1) ! mol/m2/s, tracet 0 on PE0 is mean endif c13_ocean_n=c13_ocean_n+1 endif nullify(u10m) nullify(v10m) nullify(ci) ! The rest only useful to the forward code if(enkf_is_inverse()) return ! Apply net flux scaling by ensemble member, change units, and ! coarsen to the regions. When parameters are used beyond simple ! scaling of the net flux, some piece of logic from above will ! have to be subsumed by this loop to get the parameter ! deviations. do mem=1,nmembersloc c13_ocean_mem=c13_ocean*fluxmultiplier(:,:,mem) ! parameter applied, mol/m2/s ! accumulate weekly average c13_ocean_mem=c13_ocean_mem*xmc13/1.e3 ! from mol/m2/s to kgC/m2/s call flux_to_gridbox_per_month(c13_ocean_mem) ! to kgC/area/month call coarsen_emission_2d('c13_ocean',360,180,c13_ocean_mem,c13_ocean_regxmem(:,mem),add_field) enddo end subroutine calc_emission_c13_ocean subroutine emission_apply_c13_ocean(region) ! This routine should only be called by the forward enkf code. It ! applies the ocean co2 flux by member and region to the ico2ocean ! tracer. In the inverse code, the calculated flux from each ! module is applied to the (single) c13 tracer. use dims, only: okdebug use ParTools, only: tracer_active, tracer_active use chem_param, only : ico2_oce, ic13_oce, xmc, xmc13, xmco2, xm13co2 implicit none integer, intent(in) :: region integer :: mem character(len=*), parameter :: rname = mname//'/emission_apply_c13_ocean' if(enkf_is_forward()) then if(tracer_active(ic13_oce)) then call do_add_2d(region,ic13_oce,1,c13_ocean_regxmem(region,1)%surf,xm13co2,xmc13) endif endif if(okdebug) write(*,*) 'end of emission_apply_c13_ocean' end subroutine emission_apply_c13_ocean subroutine free_emission_c13_ocean use dims, only: nregions implicit none integer :: region, n character(len=*), parameter :: rname = mname//'/free_emission_c13_ocean' if(enkf_is_forward()) then do region=1,nregions do n=1,nmembersloc deallocate(c13_ocean_regxmem(region,n)%surf) end do enddo deallocate(c13_ocean_regxmem) end if end subroutine free_emission_c13_ocean end module emission_c13_ocean