! #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" ! !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: EMISSION_DMS ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! MODULE EMISSION_DMS ! ! !USES: ! use GO, only : gol, goErr, goPr use tm5_distgrid, only : dgrid, get_distgrid, scatter, gather use partools, only : isRoot use dims, only : nregions use global_types, only : emis_data, d3_data IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: emission_dms_init ! allocate public :: emission_dms_declare ! read monthly input public :: emission_dms_apply ! distribute & add emissions to tracer array public :: emission_dms_done ! deallocate public :: getDMS ! compute DMS from ocean ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'emission_dms' type(emis_data), dimension(nregions), target :: dms_land, dms_sea ! land and sea flux type(emis_data), dimension(nregions), target :: dms_tool ! work array real, dimension(:,:), allocatable :: dms_conc ! 1x1 ocean concentrations ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - standardized routines name, new apply method ! 22 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! (1) see getDMS ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_DMS_INIT ! ! !DESCRIPTION: Allocate space needed to handle the emissions. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_DMS_INIT( status ) ! ! !USES: ! use dims, only : im, jm, iglbsfc ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - extracted from old 'declare' routine ! 22 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_DMS_Init' integer :: region integer :: imr, jmr, i1, i2, j1, j2 ! --- begin -------------------------------------- CALL GET_DISTGRID( dgrid(iglbsfc), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate( dms_conc(i1:i2, j1:j2) ) do region = 1, nregions CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr = im(region) ; jmr = jm(region) allocate(dms_land(region)%surf(i1:i2,j1:j2)) allocate(dms_sea (region)%surf(i1:i2,j1:j2)) ! work array used in coarsen/scatter if (isRoot) then allocate(dms_tool(region)%surf(imr, jmr )) else allocate(dms_tool(region)%surf(1, 1)) end if enddo status = 0 END SUBROUTINE EMISSION_DMS_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_DMS_DONE ! ! !DESCRIPTION: Deallocate space needed to handle the emissions. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_DMS_DONE( status ) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - rename old 'free_emission_dms' ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Emission_DMS_Done' integer :: region ! --- begin --------------------------------- do region = 1, nregions deallocate( dms_land (region)%surf) deallocate( dms_tool (region)%surf) deallocate( dms_sea (region)%surf) end do deallocate(dms_conc) status = 0 END SUBROUTINE EMISSION_DMS_DONE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_DMS_DECLARE ! ! !DESCRIPTION: Opens, reads and evaluates input files (per month). ! Land surface fluxes is the Spiro dataset. ! Ocean concentrations are used to retrieve fluxes by ! temperature and wind speed. !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_DMS_DECLARE( status ) ! ! !USES: ! USE MDF, ONLY : MDF_Open, MDF_HDF4, MDF_READ, MDF_Get_Var, MDF_Close use toolbox, only : coarsen_emission use dims, only : idate, sec_month, sec_year, iglbsfc, nlat180, nlon360 use chem_param, only : xms use emission_data, only : msg_emis, emis_input_dir_dms ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - updated to emission standard ! 22 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/declare_emission_dms' integer :: region, i1, i2, j1, j2, nlatsrc, nlonsrc integer, parameter :: add_field=0 integer, parameter :: amonth=2 real :: year2month integer :: hid ! work array on 1x1 resolution real, dimension(:,:), allocatable :: glb_dms_conc ! --- begin --------------------------------- year2month = sec_month/sec_year ! scale factor for yearly total to specific month nlatsrc = nlat180 ! or dgrid(iglbsfc)%jm_region, accessible thru get_distgrid(...) nlonsrc = nlon360 ! or dgrid(iglbsfc)%im_region, accessible thru get_distgrid(...) ! to deal with sea water at 1x1 if(isRoot)then allocate(glb_dms_conc(nlonsrc,nlatsrc)) else allocate(glb_dms_conc(1,1)) end if ! ! DMS emissions from landsurface ! DMS seawater concentrations. The actual flux is a function of windspeed and temperature ! and is calculated each time step. The resolution dependency needs to be checked. ! We will use the 2m winds on 1x1 resolution..... ! write(gol,'(" EMISS-INFO ------------- read DMS emissions -------------")'); call goPr if (isRoot) then ! monthly emissions from Land Spiro dataset CALL MDF_Open( trim(emis_input_dir_dms)//'/DMSland.hdf', MDF_HDF4, MDF_READ, hid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( hid, idate(2), glb_dms_conc, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Close( hid, status ) IF_NOTOK_RETURN(status=1) call msg_emis(amonth, 'SPIRO', 'vegetation/soil', 'DMS', xms, sum(glb_dms_conc)) call coarsen_emission('dms_land', nlonsrc, nlatsrc, glb_dms_conc, dms_tool, add_field, status) IF_NOTOK_RETURN(status=1) ! seawater concentrations CALL MDF_Open( trim(emis_input_dir_dms)//'/DMSconc.hdf', MDF_HDF4, MDF_READ, hid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( hid, idate(2), glb_dms_conc, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Close( hid, status ) IF_NOTOK_RETURN(status=1) end if do region = 1, nregions call scatter(dgrid(region), dms_land(region)%surf, dms_tool(region)%surf, 0, status) IF_NOTOK_RETURN(status=1) end do call scatter(dgrid(iglbsfc), dms_conc, glb_dms_conc, 0, status) IF_NOTOK_RETURN(status=1) deallocate(glb_dms_conc) ! ok status = 0 END SUBROUTINE EMISSION_DMS_DECLARE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_DMS_APPLY ! ! !DESCRIPTION: Take monthly emissions, and ! - split them vertically ! - apply time splitting factors ! - add them up (add_3d) !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_DMS_APPLY( region, status ) ! ! !USES: ! use dims, only : okdebug, itaur, nsrce, tref use dims, only : im, jm, lm use datetime, only : tau2date use emission_data, only : emission_vdist_by_sector use emission_data, only : do_add_3d use chem_param, only : idms, xmdms, xms use emission_data, only : vd_class_name_len ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - updated to use new vertical distribution method ! 28 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/emission_dms_apply' integer,dimension(6) :: idater real :: dtime, fraction integer :: imr, jmr, lmr, i1, i2, j1, j2 type(d3_data) :: emis3d character(len=vd_class_name_len) :: splittype ! --- begin ----------------------------------------- if( okdebug ) then write(gol,*) 'start of emission_dms_apply'; call goPr endif call tau2date(itaur(region),idater) dtime=float(nsrce)/(2*tref(region)) !emissions are added in two steps...XYZECCEZYX. if( okdebug ) then write(gol,*) 'emission_dms_apply in region ',region,' at date: ',idater, ' with time step:', dtime ; call goPr endif ! get a working structure for 3d emissions call get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0 ! default: no additional splitting fraction = 1.0 ! ---------------------------------------------------------------------------------------- ! distinguish here between sectors and whether they should have additional splitting ! if( ar5_sectors(lsec)%catname == 'biomassburning' ) fraction = fraction * bb_frac etc... ! ---------------------------------------------------------------------------------------- splittype = 'surface' ! ------- ! dms sea ! ------- ! vertically distribute according to sector call emission_vdist_by_sector( splittype, 'DMS', region, dms_sea(region), emis3d, status ) IF_NOTOK_RETURN(status=1;deallocate(emis3d%d3)) ! add dataset call do_add_3d( region, idms, i1, j1, emis3d%d3, xmdms, xms, status, fraction ) IF_NOTOK_RETURN(status=1) ! ------- ! dms land ! ------- emis3d%d3 = 0.0 ! vertically distribute according to sector call emission_vdist_by_sector( splittype, 'DMS', region, dms_land(region), emis3d, status ) IF_NOTOK_RETURN(status=1;deallocate(emis3d%d3)) ! add dataset call do_add_3d( region, idms, i1, j1, emis3d%d3, xmdms, xms, status, fraction ) IF_NOTOK_RETURN(status=1) if( okdebug ) then write(gol,*) 'end of emission_dms_apply'; call goPr endif deallocate( emis3d%d3 ) ! OK status = 0 END SUBROUTINE EMISSION_DMS_APPLY !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: GETDMS ! ! !DESCRIPTION: calculate the DMS flux based on high resolution ocean ! concentration and wind speeds and temperatures, ! and calculate the 'coarser' resolution fluxes from it ! ! input needed: U2 and T2 on 1x1 !\\ !\\ ! !INTERFACE: ! SUBROUTINE GETDMS( status ) ! ! !USES: ! use toolbox, only : coarsen_emission use binas, only : T0 use meteodata, only : t2m_dat, u10m_dat, v10m_dat use datetime, only : tau2date use Dims, only : itau, staggered, sec_year, sec_month, dxy11, iglbsfc, im, jm use chem_param, only : xms ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 22 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! (1) Required inputs: U2 and T2 on 1x1 ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/getdms' real, dimension(:,:), allocatable :: emi_dms, emis_glb, u_wind, v_wind, temperature real :: xwind, zschmidt, resfac, tot_flux, kw, tt, xlon, xlat, piston_dms integer :: i, j, l, region, i1, j1, i2, j2 integer, dimension(6) :: idates integer, parameter :: level1=1, add_field=0 integer :: nlatsrc, nlonsrc ! determine correction factor to account for less effective mixing ! at lower averaged windspeeds occuring at lower resolution CALL GET_DISTGRID( dgrid(iglbsfc), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) nlatsrc = jm(iglbsfc) ! or nlat180, or dgrid(iglbsfc)%jm_region, accessible thru get_distgrid(...) nlonsrc = im(iglbsfc) ! or nlon360, or dgrid(iglbsfc)%jm_region, accessible thru get_distgrid(...) if (iglbsfc /= 1) then ! assume no zoom region if region #1 is at 1x1 (see below) if (isRoot)then allocate(emis_glb (nlonsrc,nlatsrc)) else allocate(emis_glb (1,1)) end if end if allocate(emi_dms (i1:i2, j1:j2)) allocate(u_wind (i1:i2, j1:j2)) allocate(v_wind (i1:i2, j1:j2)) allocate(temperature(i1:i2, j1:j2)) call tau2date(itau+staggered,idates) !CMK at end 3-hourly interval temperature = t2m_dat(iglbsfc)%data(:,:,1) ! convert to celsius temperature = temperature - T0 u_wind = u10m_dat(iglbsfc)%data(:,:,1) v_wind = v10m_dat(iglbsfc)%data(:,:,1) resfac=1. ! THE RESOLUTION DEPENDENCY WITH PREVIOUS TM3 RUNS WAS AS FOLLOWS ! NEEDS TO BE CHECKED FD !FD if (jm.eq.24) resfac=1.32 !FD if (jm.eq.48) resfac=1.065 emi_dms=0. do j=j1,j2 do i=i1,i2 ! ! ocean surface temperature is approximated by air surface temperature ! we assume that if temperature<-20. ! sea ice prevents DMS emissions (sea-ice data from ECMWF would be usefull) ! max ocean temperature is 28.0 C ! minumum ocean temperature 0. C !fd tt=min(28.,ts(i,j)-273.15) ! the original parameterisation needs the seawater temperature ! as a proxy we use the air temperature, but restrict it to 28 C. ! surface or 2m temperature of air is probably the best proxy. tt=min(28.0,temperature(i,j)) if (dms_conc(i,j).gt.0..and.tt.gt.-20.) then zschmidt=max(0.0,3652.047271-246.99*TT+8.536397*TT*TT-0.124397*TT*TT*TT) xwind=sqrt(u_wind(i,j)**2+v_wind(i,j)**2) ! ! Liss and Merlivat windspeed dependency of transfer velocity kw ! including correction for Schmidt number ! windspeed at middle of gridbox is used should be 10 m windspeed [m/s] if (xwind.le.3.6) then kw=0.17*xwind piston_dms=kw*(zschmidt/600.)**(-0.666) else !xwind>3.6) if (xwind.le.13.0) then kw=2.85*xwind-9.65 piston_dms=kw*(zschmidt/600.)**(-0.5) else !xwind>13.0 kw=5.9*xwind-49.3 piston_dms=kw*(zschmidt/600.)**(-0.5) endif !x<13.0 endif ! ! cm/hr=> m/s * nmol/l=>mol/m3 * m2 * g S/mol=>g S/s ! emi_dms(i,j)=resfac*piston_dms/360000.*dms_conc(i,j)*1e-6*dxy11(j)*xms endif !(dms_conc gt 0) enddo enddo !----------------------------------------------------------- ! Fill target array !----------------------------------------------------------- ! Test for optimization: if region #1 is at 1x1, assume no zoom region and skip call to coarsen and mpi comm if (iglbsfc == 1) then dms_sea(1)%surf = emi_dms*sec_month*1e-3 !expected in kgS/month else !----------------------------------------------------------- ! gather on 1x1, coarsen to other grid on root, then scatter !----------------------------------------------------------- ! FIXME-PERF : Need a coarsen that works independtly on each proc call gather(dgrid(iglbsfc), emi_dms, emis_glb, 0, status) IF_NOTOK_RETURN(status=1) if (isRoot) then ! FIXME : if really needed, write into budget - here, it gives too much output !write(gol,*)'the yearly sum of GETDMS would be ',sum(emis_glb) *1e-12*sec_year,' Tg S/a'; call goPr ! assume equal emissions over season call coarsen_emission('emi_dms', nlonsrc, nlatsrc, emis_glb, dms_tool, add_field, status) IF_NOTOK_RETURN(status=1) end if do region = 1, nregions call scatter(dgrid(region), dms_sea(region)%surf, dms_tool(region)%surf, 0, status) IF_NOTOK_RETURN(status=1) dms_sea(region)%surf = dms_sea(region)%surf*sec_month*1e-3 !expected in kgS/month enddo deallocate(emis_glb) end if !----------------------------------------------------------- deallocate(emi_dms) deallocate(temperature) deallocate(u_wind) deallocate(v_wind) END SUBROUTINE GETDMS !EOC END MODULE EMISSION_DMS