diff --git a/atmos_shared/atmos_cmip_diag.F90 b/atmos_shared/atmos_cmip_diag.F90 new file mode 100644 index 00000000..4b0e91cc --- /dev/null +++ b/atmos_shared/atmos_cmip_diag.F90 @@ -0,0 +1,763 @@ +module atmos_cmip_diag_mod + +!---------------------------------------------------------------------- +! Module for registering and sending (writing) 3D CMIP diagnostic +! data on model levels and pressure levels. New vertical axes are +! defined from lowest to uppermost level (opposite the model). +! Prefined pressure axes corresponding to CMIP axes are used. +! The vertical axis used is specified via the 'module_name' field +! in the diag_table. +!---------------------------------------------------------------------- + +use mpp_mod, only: input_nml_file +use fms_mod, only: check_nml_error, & + stdlog, mpp_pe, mpp_root_pe, & + write_version_number, & + error_mesg, FATAL, WARNING, NOTE, & + lowercase, string +use time_manager_mod, only: time_type +use diag_manager_mod, only: diag_axis_init, register_diag_field, & + send_data, get_diag_field_id, & + register_static_field, & + diag_axis_add_attribute, & + diag_field_add_attribute, & + DIAG_FIELD_NOT_FOUND +use diag_data_mod, only: CMOR_MISSING_VALUE + +!---------------------------------------------------------------------- + +implicit none +private + +!---------------------------------------------------------------------- + +public :: atmos_cmip_diag_init, atmos_cmip_diag_end, & + register_cmip_diag_field_2d, & + register_cmip_diag_field_3d, & + send_cmip_data_3d, & + query_cmip_diag_id + +!---------------------------------------------------------------------- + +! vertical pressure grids +! plev = same as plev17, unless use_extra_levels = .true. +! plev19 = Table Amon = standard 17 levels + 5, 1 hPa +! plev8 = daily data +! plev7h = HighResMIP (6hr time mean, 3hr synoptic) +! plev3 = used in CMIP5 for 6hrPlev +! plev23 = plev19 + (7,3,2,0.4 hPa) +! plev4 = used in 6hrPlev + +real, dimension(39) :: plev39 = & + (/ 100000., 92500., 85000., 70000., 60000., 50000., & + 40000., 30000., 25000., 20000., 17000., 15000., & + 13000., 11500., 10000., 9000., 8000., 7000., & + 5000., 3000., 2000., 1500., 1000., 700., & + 500., 300., 200., 150., 100., 70., & + 50., 40., 30., 20., 15., 10., & + 7., 5., 3. /) +real, dimension(26) :: plev26 = & + (/ 100000., 92500., 85000., 70000., 60000., 50000., & + 40000., 30000., 25000., 20000., 17000., 15000., & + 13000., 11500., 10000., 9000., 8000., 7000., & + 5000., 3000., 2000., 1500., 1000., 700., & + 500., 300. /) +real, dimension(23) :: plev23 = & + (/ 100000., 92500., 85000., 70000., 60000., 50000., & + 40000., 30000., 25000., 20000., 15000., 10000., & + 7000., 5000., 3000., 2000., 1000., 700., & + 500., 300., 200., 100., 40. /) +real, dimension(19) :: plev19 = & + (/ 100000., 92500., 85000., 70000., 60000., 50000., & + 40000., 30000., 25000., 20000., 15000., 10000., & + 7000., 5000., 3000., 2000., 1000., 500., & + 100. /) +real, dimension(8) :: plev8 = & + (/ 100000., 85000., 70000., 50000., & + 25000., 10000., 5000., 1000. /) +real, dimension(7) :: plev7h = & + (/ 92500., 85000., 70000., 60000., 50000., 25000., 5000. /) +real, dimension(4) :: plev4 = & + (/ 92500., 85000., 50000., 25000. /) +real, dimension(3) :: plev3 = & + (/ 85000., 50000., 25000. /) + +!----------------------------------------------------------------------- +!--- namelist --- + +logical :: use_extra_levels = .true. ! use more than the standard + ! 17 pressure levels when possible + +logical :: flip_cmip_levels = .true. ! flip vertical model level output + ! from bottom(surface) to top. + +logical :: output_modeling_realm = .false. ! add modeling_realm attribute + ! to all variables + +logical :: error_when_phalf_missing = .true. ! when pressure level interp is requested + ! and log(phalf) is not supplied a fatal + ! error will result, set to false to get + ! the previous behavior (warsaw_201803). + +character(len=64) :: modeling_realm_default = 'atmos' ! default modeling_realm attribute + ! can be overriden in + ! register_cmip_diag + +integer :: verbose = 1 ! verbose level = 0,1,2 + +namelist /atmos_cmip_diag_nml/ use_extra_levels, flip_cmip_levels, & + output_modeling_realm, modeling_realm_default, & + error_when_phalf_missing, verbose + +!----------------------------------------------------------------------- + +integer, parameter :: MAXPLEVS = 9 ! max plev sets +integer, dimension(MAXPLEVS) :: num_pres_levs +real, dimension(MAXPLEVS,50) :: pressure_levels ! max 50 levels per set + +character(len=16) :: mod_name = 'cmip' + +! cmip vertical axis names +! index -1 = 'levhalf' half model levels +! index 0 = 'lev' full model levels +! index >0 = 'plev*' pressure levels +character(len=128), dimension(-1:MAXPLEVS) :: cmip_axis_names +integer, dimension(3,-1:MAXPLEVS) :: cmip_axis_data + +integer :: area_id + +!---------------------------------------------------------------------- + +!--- store field id for all possible axes +! index 0 = on model level (either full or half) +! index >0 = on pressure levels +public :: cmip_diag_id_type +type cmip_diag_id_type + integer, dimension(0:MAXPLEVS) :: field_id = 0 +end type cmip_diag_id_type + +!---------------------------------------------------------------------- + +character(len=128) :: version = '$Id$' +character(len=128) :: tagname = '$Name$' + +logical :: module_is_initialized=.false. + +CONTAINS + +!####################################################################### + +subroutine atmos_cmip_diag_init ( ak, bk, ptop, axes, Time ) +real, intent(in), dimension(:) :: ak, bk ! ap,b at model layer interfaces +real, intent(in) :: ptop ! pressure at top level +integer, intent(in) :: axes(2) ! x/y axes identifiers +type(time_type), intent(in) :: Time + +!----------------------------------------------------------------------- +! local data + +integer :: axes3d(3), k, kk, ind, np, id_plev, num_std_plevs +integer :: nlev +integer :: iunit, ierr, io +logical :: used +character(len=16) :: axis_name +integer :: flip + +real :: p0 +real, dimension(size(ak,1)-1) :: ap, b, lev +real, dimension(2,size(ak,1)-1) :: ap_bnds, b_bnds, lev_bnds +real, dimension(size(ak,1)) :: levhalf + +integer :: id_lev, id_levhalf, id_nv, id_ap, id_b, & + id_ap_bnds, id_b_bnds, id_lev_bnds + +!----------------------------------------------------------------------- + + if (module_is_initialized) then + call error_mesg ('atmos_cmip_diag_mod', & + 'module has already been initialized', WARNING) + return + endif + +!----------------------------------------------------------------------- +!----- read namelist ----- + read (input_nml_file, nml=atmos_cmip_diag_nml, iostat=io) + ierr = check_nml_error (io, 'atmos_cmip_diag_nml') + +!----- write version and namelist to log file ----- + + iunit = stdlog() + call write_version_number ( version, tagname ) + if (mpp_pe() == mpp_root_pe()) write (iunit, nml=atmos_cmip_diag_nml) + + +!----------------------------------------------------------------------- +! axis and area identifiers + axes3d(1:2) = axes(1:2) + area_id = get_diag_field_id ('dynamics', 'area') + if (area_id .eq. DIAG_FIELD_NOT_FOUND) call error_mesg & + ('atmos_cmip_diag_init', 'diagnostic field "dynamics", '// & + '"area" is not in the diag_table', NOTE) + +!----------------------------------------------------------------------- +! determine the maximum number of standard pressure levels +! first get the pressure (based on ps=1000hPa) at top model level + + if (use_extra_levels) then + do k = 23, 1, -1 + if (plev23(k) .gt. ptop) then + num_std_plevs = k + exit + endif + enddo + else + num_std_plevs = 17 ! standard ncep 17 levels + endif + +!----------------------------------------------------------------------- +! vertical coordinate variables +! cmip levels are defined from the surface to top (flip model values) + + flip = 1 + if (flip_cmip_levels) flip = -1 + nlev = size(ak,1)-1 + + p0 = 100000. + do k = 1, nlev + if (flip_cmip_levels) then + kk = nlev - k + 2 + else + kk = k + end if + ap_bnds(1,k) = ak(kk) + ap_bnds(2,k) = ak(kk+flip) + b_bnds(1,k) = bk(kk) + b_bnds(2,k) = bk(kk+flip) + ap(k) = (ap_bnds(1,k)+ap_bnds(2,k))*0.5 + b(k) = (b_bnds(1,k)+b_bnds(2,k))*0.5 + enddo + lev = ap/p0 + b ! definition for CMIP purposes + lev_bnds = ap_bnds/p0 + b_bnds + levhalf(1:nlev) = ap_bnds(1,:)/p0 + b_bnds(1,:) ! definition at half levels for CMIP purposes + levhalf(nlev+1) = ap_bnds(2,nlev)/p0 + b_bnds(2,nlev) + + !---- register new axes ---- + + ! at full levels (with bounds attribute) + id_lev = diag_axis_init('lev', lev, '1.0', 'Z', & + 'hybrid sigma pressure coordinate', & + direction=-1, set_name='cmip', req='lev_bnds') + call diag_axis_add_attribute (id_lev, 'formula', 'p(n,k,j,i) = ap(k) + b(k)*ps(n,j,i)') + call diag_axis_add_attribute (id_lev, 'formula_terms', 'ap: ap b: b ps: ps') + call diag_axis_add_attribute (id_lev, 'bounds', 'lev_bnds') + call diag_axis_add_attribute (id_lev, 'standard_name', & + 'atmosphere_hybrid_sigma_pressure_coordinate') + + ! at half levels (bounds unknown at top and bottom) + id_levhalf = diag_axis_init('levhalf', levhalf, '1.0', 'Z', & + 'hybrid sigma pressure coordinate', & + direction=-1, set_name='cmip') + call diag_axis_add_attribute (id_levhalf, 'standard_name', & + 'atmosphere_hybrid_sigma_pressure_coordinate') + call diag_axis_add_attribute ( id_levhalf, 'formula', 'p(n,k+1/2,j,i) = ap(k+1/2) + b(k+1/2)*ps(n,j,i)') + call diag_axis_add_attribute ( id_levhalf, 'formula_terms', 'ap: ap_bnds b: b_bnds ps: ps') + + ! vertex number for bounds dimension + id_nv = diag_axis_init('nv', (/1.,2./), 'none', 'N', 'vertex number', set_name='nv') + + ! register new static variables + + id_ap = register_static_field (mod_name, 'ap', (/id_lev/), & + 'vertical coordinate formula term: ap(k)', 'Pa') + + id_b = register_static_field (mod_name, 'b', (/id_lev/), & + 'vertical coordinate formula term: b(k)', '1.0') + + id_ap_bnds = register_static_field (mod_name, 'ap_bnds', (/id_nv,id_lev/), & + 'vertical coordinate formula term: ap(k+1/2)', 'Pa') + + id_b_bnds = register_static_field (mod_name, 'b_bnds', (/id_nv,id_lev/), & + 'vertical coordinate formula term: b(k+1/2)', '1.0') + + id_lev_bnds = register_static_field (mod_name, 'lev_bnds', (/id_nv,id_lev/), & + 'hybrid sigma pressure coordinate', '1.0', & + standard_name='atmosphere_hybrid_sigma_pressure_coordinate') + if (id_lev_bnds > 0) then + call diag_field_add_attribute ( id_lev_bnds, 'formula', 'p(n,k+1/2,j,i) = ap(k+1/2) + b(k+1/2)*ps(n,j,i)') + call diag_field_add_attribute ( id_lev_bnds, 'formula_terms', 'ap: ap_bnds b: b_bnds ps: ps') + endif + + ! save static data + if (id_ap > 0) used = send_data ( id_ap, ap, Time ) + if (id_b > 0) used = send_data ( id_b , b , Time ) + if (id_ap_bnds > 0) used = send_data ( id_ap_bnds, ap_bnds, Time ) + if (id_b_bnds > 0) used = send_data ( id_b_bnds, b_bnds, Time ) + if (id_lev_bnds > 0) used = send_data ( id_lev_bnds, lev_bnds, Time ) + + axes3d(3) = id_lev + cmip_axis_names(0) = 'lev' !mod_name + cmip_axis_data(:,0) = axes3d + + axes3d(3) = id_levhalf + cmip_axis_names(-1) = 'levhalf' !mod_name + cmip_axis_data(:,-1) = axes3d +!----------------------------------------------------------------------- +! loop through all possible pressure level sets +! initialize the pressure axis +! define new 3d grid +! define all 3d state variable on this 3d grid + + do ind = 1, MAXPLEVS + if (ind .eq. 1) then + np = num_std_plevs + pressure_levels(ind,1:np) = plev23(1:np) + axis_name = 'plev_std' + else if (ind .eq. 2) then + np = size(plev19,1) + pressure_levels(ind,1:np) = plev19 + axis_name = 'plev19' + else if (ind .eq. 3) then + np = size(plev8,1) + pressure_levels(ind,1:np) = plev8 + axis_name = 'plev8' + else if (ind .eq. 4) then + np = size(plev3,1) + pressure_levels(ind,1:np) = plev3 + axis_name = 'plev3' + else if (ind .eq. 5) then + np = size(plev7h,1) + pressure_levels(ind,1:np) = plev7h + axis_name = 'plev7h' + else if (ind .eq. 6) then + np = size(plev23,1) + pressure_levels(ind,1:np) = plev23 + axis_name = 'plev23' + else if (ind .eq. 7) then + np = size(plev4,1) + pressure_levels(ind,1:np) = plev4 + axis_name = 'plev4' + else if (ind .eq. 8) then + np = size(plev26,1) + pressure_levels(ind,1:np) = plev26 + axis_name = 'plev26' + else if (ind .eq. 9) then + np = size(plev39,1) + pressure_levels(ind,1:np) = plev39 + axis_name = 'plev39' + endif + + num_pres_levs(ind) = np + id_plev = diag_axis_init(axis_name, pressure_levels(ind,1:np), & + 'Pa', 'z', 'pressure', direction=-1, set_name="dynamics") + + axes3d(3) = id_plev + cmip_axis_names(ind) = trim(axis_name) !trim(mod_name)//'_'//trim(axis_name) + cmip_axis_data(:,ind) = axes3d + + enddo + + if (verbose > 0) then + call error_mesg('atmos_cmip_diag_mod', & + 'cmip_axis_names(-1) = "'//trim(cmip_axis_names(-1))//'"',NOTE) + do ind = 0, MAXPLEVS + call error_mesg('atmos_cmip_diag_mod', & + 'cmip_axis_names('//trim(string(ind))//') = "'//trim(cmip_axis_names(ind))//'"',NOTE) + enddo + endif + +!--- done --- + module_is_initialized=.true. + +!----------------------------------------------------------------------- + +end subroutine atmos_cmip_diag_init + +!####################################################################### + +logical function query_cmip_diag_id (cmip_id, pres) +type(cmip_diag_id_type), intent(in) :: cmip_id +logical, optional, intent(in) :: pres +integer :: is + + is = 0 + if (present(pres)) then + if (pres) is = 1 + end if + + query_cmip_diag_id = count(cmip_id%field_id(is:) > 0) .gt. 0 + +end function query_cmip_diag_id + +!####################################################################### + +integer function register_cmip_diag_field_2d (module_name, field_name, & + Time_init, long_name, units, standard_name, & + missing_value, interp_method, mask_variant, realm) + + character(len=*), intent(in) :: module_name, field_name + type(time_type), intent(in) :: Time_init + character(len=*), intent(in), optional :: long_name, units, standard_name + real, intent(in), optional :: missing_value + character(len=*), intent(in), optional :: interp_method, realm + logical , intent(in), optional :: mask_variant + + real :: mvalue + character(len=64) :: modeling_realm +!----------------------------------------------------------------------- + mvalue = CMOR_MISSING_VALUE; if (present(missing_value)) mvalue = missing_value + modeling_realm = modeling_realm_default + if (present(realm)) modeling_realm = realm + + if (output_modeling_realm) then + register_cmip_diag_field_2d = register_diag_field (module_name, field_name, & + cmip_axis_data(1:2,0), Time_init, long_name=long_name, & + units=units, standard_name=standard_name, area=area_id, & + mask_variant=mask_variant, missing_value=mvalue, & + interp_method=interp_method, realm=modeling_realm ) + else + register_cmip_diag_field_2d = register_diag_field (module_name, field_name, & + cmip_axis_data(1:2,0), Time_init, long_name=long_name, & + units=units, standard_name=standard_name, area=area_id, & + mask_variant=mask_variant, missing_value=mvalue, & + interp_method=interp_method) + endif + +!----------------------------------------------------------------------- + +end function register_cmip_diag_field_2d + +!####################################################################### + +function register_cmip_diag_field_3d (module_name, field_name, & + Time_init, long_name, units, standard_name, & + axis, missing_value, interp_method, mask_variant, & + realm) + + character(len=*), intent(in) :: module_name, field_name + type(time_type), intent(in) :: Time_init + character(len=*), intent(in), optional :: long_name, units, standard_name + real, intent(in), optional :: missing_value + character(len=*), intent(in), optional :: axis ! 'full' or 'half' levels + character(len=*), intent(in), optional :: interp_method ! for fregrid + logical , intent(in), optional :: mask_variant + character(len=*), intent(in), optional :: realm ! modeling realm + + type(cmip_diag_id_type) :: register_cmip_diag_field_3d + integer :: ind, indx, kount + real :: mvalue + character(len=128) :: module_name_table + character(len=4) :: vert_axis + character(len=64) :: modeling_realm +!----------------------------------------------------------------------- + + mvalue = CMOR_MISSING_VALUE; if (present(missing_value)) mvalue = missing_value + vert_axis = 'full'; if (present(axis)) vert_axis = lowercase(trim(axis)) + + modeling_realm = modeling_realm_default + if (present(realm)) modeling_realm = realm + + register_cmip_diag_field_3d%field_id = 0 + + ! loop thru all axes + do ind = 0, MAXPLEVS + indx = ind + if (ind .eq. 0 .and. vert_axis .eq. 'half') indx = -1 + + module_name_table = trim(module_name) + if (ind .gt. 0) then + module_name_table = trim(module_name_table)//'_'//trim(cmip_axis_names(ind)) + end if + + ! only register fields that are in the diag_table + if ( get_diag_field_id(module_name_table, field_name) .ne. DIAG_FIELD_NOT_FOUND ) then + + if (output_modeling_realm) then + register_cmip_diag_field_3d%field_id(ind) = register_diag_field(module_name_table, field_name, & + cmip_axis_data(:,indx), Time_init, long_name=long_name, units=units, & + standard_name=standard_name, area=area_id, mask_variant=mask_variant, & + missing_value=mvalue, interp_method=interp_method, realm=modeling_realm) + else + register_cmip_diag_field_3d%field_id(ind) = register_diag_field(module_name_table, field_name, & + cmip_axis_data(:,indx), Time_init, long_name=long_name, units=units, & + standard_name=standard_name, area=area_id, mask_variant=mask_variant, & + missing_value=mvalue, interp_method=interp_method) + endif + + if (verbose > 0) call error_mesg('atmos_cmip_diag_mod', & + 'register cmip diag: module='//trim(module_name_table)//', field='//trim(field_name)// & + ', field_id='//trim(string(register_cmip_diag_field_3d%field_id(ind))),NOTE) + + else if (verbose > 1) then + ! for debugging purposes + call error_mesg('atmos_cmip_diag_mod','NOT registering cmip diag: module='// & + trim(module_name_table)//', field='//trim(field_name),NOTE) + endif + enddo + + if (verbose > 1) then + kount = count(register_cmip_diag_field_3d%field_id > 0) + if (query_cmip_diag_id(register_cmip_diag_field_3d)) then + call error_mesg('atmos_cmip_diag_mod','query_cmip_diag_id=TRUE, module='// & + trim(module_name_table)//', field='//trim(field_name)//', kount='//trim(string(kount)),NOTE) + else + call error_mesg('atmos_cmip_diag_mod','query_cmip_diag_id=FALSE, module='// & + trim(module_name_table)//', field='//trim(field_name)//', kount='//trim(string(kount)),NOTE) + endif + endif + +!----------------------------------------------------------------------- + +end function register_cmip_diag_field_3d + +!####################################################################### + +logical function send_cmip_data_3d (cmip_id, field, Time, is_in, js_in, ks_in, phalf, mask, rmask, opt, ext) + + type(cmip_diag_id_type), intent(in) :: cmip_id + real, dimension(:,:,:), intent(in) :: field + type(time_type), intent(in), optional :: Time + integer, intent(in), optional :: is_in, js_in, ks_in + real, dimension(:,:,:), intent(in), optional :: phalf, rmask + logical, dimension(:,:,:), intent(in), optional :: mask + integer, intent(in), optional :: opt ! if opt /= 0 then phalf(i,k,j) + logical, intent(in), optional :: ext + + integer :: ind, id, np, ke + real, allocatable :: pdat(:,:,:) + +!----------------------------------------------------------------------- + + if (.not.module_is_initialized) call error_mesg ('atmos_cmip_diag_mod', & + 'module has not been initialized', FATAL) + + if (present(ks_in)) then + if (ks_in .ne. 1) call error_mesg ('atmos_cmip_diag_mod', & + 'subroutine send_cmip_data_3d does not support optional arg "ks_in"', FATAL) + endif + + if (present(rmask) .and. present(mask)) call error_mesg('atmos_cmip_diag_mod', & + 'rmask and mask can not both be present',FATAL) + + send_cmip_data_3d = .false. + + ! loop thru all axes + do ind = 0, MAXPLEVS + + if (cmip_id%field_id(ind) > 0) then + id = cmip_id%field_id(ind) + + ! pressure level interpolation if "phalf" is present + + ! pressure level interpolation when ind > 0 + if (ind > 0) then + if (.not.present(phalf)) then + if (error_when_phalf_missing) then + call error_mesg('atmos_cmip_diag_mod', & + 'log(phalf) must be present for pressure level interpolation',FATAL) + else + cycle ! silently skip? + endif + endif + if (present(rmask) .or. present(mask)) call error_mesg('atmos_cmip_diag_mod', & + 'rmask or mask not allowed with pressure interpolation',FATAL) + np = num_pres_levs(ind) + allocate(pdat(size(field,1),size(field,2),np)) + call interpolate_vertical_3d (pressure_levels(ind,1:np), phalf, field, pdat, opt=opt, ext=ext) + send_cmip_data_3d = send_data(id, pdat, Time, is_in=is_in, js_in=js_in, ks_in=ks_in) + deallocate(pdat) + + else + ! save data on model levels (flip data) + if (flip_cmip_levels) then + ke = size(field,3) + if (.not.present(mask) .and. .not.present(rmask)) then + send_cmip_data_3d = send_data(id, field(:,:,ke:1:-1), Time, & + is_in=is_in, js_in=js_in, ks_in=ks_in) + else if (present(mask) .and. .not.present(rmask)) then + send_cmip_data_3d = send_data(id, field(:,:,ke:1:-1), Time, & + is_in=is_in, js_in=js_in, ks_in=ks_in, mask=mask(:,:,ke:1:-1)) + else if (.not.present(mask) .and. present(rmask)) then + send_cmip_data_3d = send_data(id, field(:,:,ke:1:-1), Time, & + is_in=is_in, js_in=js_in, ks_in=ks_in, rmask=rmask(:,:,ke:1:-1)) + endif + else + send_cmip_data_3d = send_data(id, field(:,:,:), Time, & + is_in=is_in, js_in=js_in, ks_in=ks_in, mask=mask, rmask=rmask) + endif + endif + else + send_cmip_data_3d = .false. + endif + enddo + +!----------------------------------------------------------------------- + +end function send_cmip_data_3d + +!####################################################################### + +subroutine atmos_cmip_diag_end + +! do nothing, no way to unregister diag fields + +end subroutine atmos_cmip_diag_end + +!####################################################################### + +subroutine dealloc_cmip_diag_id_type (cmip_id) +class(cmip_diag_id_type), intent(inout) :: cmip_id + + !deallocate(cmip_id%field_id) + +end subroutine dealloc_cmip_diag_id_type + +!####################################################################### +! wrapper for different vertical interpolation routines +! opt = 0 for standard indexing of peln(i,j,k) -- this is the default +! opt /= 0 for FV-core indexing of peln(i,k,j) +! ext = flag to extrapolate data below (and above) input data (default: false) + +subroutine interpolate_vertical_3d (plev, peln, a, ap, opt, ext) + real, intent(in), dimension(:) :: plev ! target p-levels + real, intent(in), dimension(:,:,:) :: peln ! log(phalf), model half levels + real, intent(in), dimension(:,:,:) :: a ! input data + real, intent(out), dimension(:,:,:) :: ap ! output data on p-levels + integer, intent(in), optional :: opt ! peln indexing + logical, intent(in), optional :: ext ! extrapolate? + + integer :: im, jm, km, kp + integer :: iopt + + iopt = 0; if (present(opt)) iopt = opt + im = size(a,1) + jm = size(a,2) + km = size(a,3) + kp = size(plev,1) + + if (iopt .eq. 0) then + if (size(peln,2).eq.jm .and. size(peln,3).eq.km+1) then + call interpolate_vertical (im, jm, km, kp, plev, peln, a, ap) ! peln(im,jm,km+1) + !else if (size(peln,2).eq.jm .and. size(peln,3).eq.km) then + ! call interpolate_vertical_half (im, jm, km, kp, plev, peln, a, ap, ext) ! peln(im,jm,km) + else + call error_mesg('atmos_cmip_diag_mod','invalid indexing option and/or array sizes',FATAL) + endif + + else + if (size(peln,3).eq.jm .and. size(peln,2).eq.km+1) then + call interpolate_vertical_fv (im, jm, km, kp, plev, peln, a, ap) ! peln(im,km+1,jm) + else if (size(peln,3).eq.jm .and. size(peln,2).eq.km) then + call interpolate_vertical_half_fv (im, jm, km, kp, plev, peln, a, ap, ext) ! peln(im,km,jm) + else + call error_mesg('atmos_cmip_diag_mod','invalid indexing option and/or array sizes',FATAL) + endif + + endif + +end subroutine interpolate_vertical_3d + +!####################################################################### +! a (im, jm, km ) <-- input data on FULL model levels +! peln (im, jm, km+1) <-- standard indexing (i,j,k) +! km = number of FULL levels + +subroutine interpolate_vertical (im, jm, km, np, plev, peln, a, ap, ext) + integer, intent(in) :: im, jm, km, np + real, intent(in), dimension(np) :: plev ! target p-levels + real, intent(in), dimension(im,jm,km+1) :: peln ! log(phaf), model half levels + real, intent(in), dimension(im,jm,km) :: a ! input data on model levels + real, intent(out), dimension(im,jm,np) :: ap ! output data on p-levels + logical, intent(in), optional :: ext + + real, dimension(km,im,jm) :: pm + real :: logp + integer :: i, j, k, kp + logical :: extrap + + extrap = .false.; if (present(ext)) extrap = ext + + do k = 1, km + do j = 1, jm + do i = 1, im + pm(k,i,j) = 0.5*(peln(i,j,k)+peln(i,j,k+1)) + enddo + enddo + enddo + + include "atmos_cmip_interp.inc" + + +end subroutine interpolate_vertical + +!####################################################################### +! a (im, jm, km) <-- input data on FULL model levels +! peln (im, km+1, jm) <-- FV core indexing +! km = number of FULL levels + +subroutine interpolate_vertical_fv (im, jm, km, np, plev, peln, a, ap, ext) + integer, intent(in) :: im, jm, km, np + real, intent(in), dimension(np) :: plev ! target p-levels + real, intent(in), dimension(im,km+1,jm) :: peln ! log(phaf), model half levels + real, intent(in), dimension(im,jm,km) :: a ! input data on model levels + real, intent(out), dimension(im,jm,np) :: ap ! output data on p-levels + logical, intent(in), optional :: ext + + real, dimension(km,im,jm) :: pm + real :: logp + integer :: i, j, k, kp + logical :: extrap + + extrap = .false.; if (present(ext)) extrap = ext + + + do j = 1, jm + do k = 1, km + do i = 1, im + pm(k,i,j) = 0.5*(peln(i,k,j)+peln(i,k+1,j)) + enddo + enddo + enddo + + include "atmos_cmip_interp.inc" + + +end subroutine interpolate_vertical_fv + +!####################################################################### +! a (im, jm, km) <-- input data on HALF model levels +! peln (im, km, jm) <-- FV core indexing +! km = number of HALF levels + +subroutine interpolate_vertical_half_fv (im, jm, km, np, plev, peln, a, ap, ext) + + integer, intent(in) :: im, jm, km, np + real, intent(in), dimension(np) :: plev ! target p-levels + real, intent(in), dimension(im,km,jm) :: peln ! log(phaf), model half levels + real, intent(in), dimension(im,jm,km) :: a ! input data on model HALF levels + real, intent(out), dimension(im,jm,np) :: ap ! output data on p-levels + logical, intent(in), optional :: ext + + real, dimension(km,im,jm) :: pm + real :: logp + integer :: i, j, k, kp + logical :: extrap + + extrap = .false.; if (present(ext)) extrap = ext + + do j = 1, jm + do k = 1, km + do i = 1, im + pm(k,i,j) = peln(i,k,j) + enddo + enddo + enddo + + include "atmos_cmip_interp.inc" + + +end subroutine interpolate_vertical_half_fv + +!####################################################################### + +end module atmos_cmip_diag_mod + diff --git a/atmos_shared/atmos_cmip_interp.inc b/atmos_shared/atmos_cmip_interp.inc new file mode 100644 index 00000000..afc699c6 --- /dev/null +++ b/atmos_shared/atmos_cmip_interp.inc @@ -0,0 +1,44 @@ + + if (extrap) then + do kp = 1, np + logp = log(plev(kp)) + do j = 1, jm + do i = 1, im + if (logp <= pm(1,i,j)) then + ap(i,j,kp) = a(i,j,1) + (a(i,j,3)-a(i,j,1)) * (logp-pm(1,i,j))/(pm(3,i,j)-pm(1,i,j)) + else if (logp >= pm(km,i,j)) then + ap(i,j,kp) = a(i,j,km) + (a(i,j,km)-a(i,j,km-2)) * (logp-pm(km,i,j))/(pm(km,i,j)-pm(km-2,i,j)) + else + do k = 1, km-1 + if (logp >= pm(k,i,j) .and. logp <= pm(k+1,i,j)) then + ap(i,j,kp) = a(i,j,k) + (a(i,j,k+1)-a(i,j,k)) * (logp-pm(k,i,j))/(pm(k+1,i,j)-pm(k,i,j)) + goto 1000 + endif + enddo + endif +1000 continue + enddo + enddo + enddo + else + do kp = 1, np + logp = log(plev(kp)) + do j = 1, jm + do i = 1, im + if (logp <= pm(1,i,j)) then + ap(i,j,kp) = a(i,j,1) + else if (logp >= pm(km,i,j)) then + ap(i,j,kp) = a(i,j,km) + else + do k = 1, km-1 + if (logp >= pm(k,i,j) .and. logp <= pm(k+1,i,j)) then + ap(i,j,kp) = a(i,j,k) + (a(i,j,k+1)-a(i,j,k)) * (logp-pm(k,i,j))/(pm(k+1,i,j)-pm(k,i,j)) + goto 2000 + endif + enddo + endif +2000 continue + enddo + enddo + enddo + endif diff --git a/atmos_shared/atmos_global_diag.F90 b/atmos_shared/atmos_global_diag.F90 new file mode 100644 index 00000000..fbe7fd18 --- /dev/null +++ b/atmos_shared/atmos_global_diag.F90 @@ -0,0 +1,430 @@ +module atmos_global_diag_mod + +!---------------------------------------------------------------------- +! Module for computed globally averaged atmospheric quantities +! and then registering and sending this data to the diag_manager. +!---------------------------------------------------------------------- + +use mpp_mod, only: input_nml_file +use mpp_domains_mod, only: domain2d, mpp_global_sum, BITWISE_EFP_SUM, & + null_domain2d, operator(.eq.) +use fms_mod, only: check_nml_error, & + stdlog, mpp_pe, mpp_root_pe, & + write_version_number, & + error_mesg, FATAL, WARNING +use time_manager_mod, only: time_type +use diag_manager_mod, only: register_diag_field, send_data, & + DIAG_FIELD_NOT_FOUND +!use diag_data_mod, only: null_axis_id, & +use diag_data_mod, only: CMOR_MISSING_VALUE +use diag_axis_mod, only: get_domain2d + +!----------------------------------------------------------------------- + +implicit none +private + +!----------------------------------------------------------------------- +! public interfaces + +public :: atmos_global_diag_init, & + register_global_diag_field, & + get_global_diag_field_id, & + buffer_global_diag, & + send_global_diag, & + atmos_global_diag_end + +!----------------------------------------------------------------------- + +interface send_global_diag + module procedure send_global_diag_data + module procedure send_global_diag_buffer +end interface + +!----------------------------------------------------------------------- +! private data structure + +type atmos_global_diag_type + character(len=128) :: field_name + integer :: field_id + real, pointer :: buffer(:,:) + type(time_type) :: Time + logical :: use_buffer + logical :: use_masking +end type atmos_global_diag_type + +!----------------------------------------------------------------------- +! private module data + +real, allocatable, dimension(:,:) :: area_g +real :: area_g_sum +integer :: id, jd +!integer :: null_axis(1) +real :: missing_value = CMOR_MISSING_VALUE + +character(len=12) :: mod_name = 'atmos_global' +type(domain2d), save :: Domain2 + +integer :: num_fields = 0 +type(atmos_global_diag_type), allocatable :: fields(:) + +character(len=128) :: version = '$Id$' +character(len=128) :: tagname = '$Name$' + +logical :: module_is_initialized=.false. + +!----------------------------------------------------------------------- +! namelist + +integer :: max_fields = 100 + +namelist /atmos_global_diag_nml/ max_fields + + +CONTAINS + +!####################################################################### + +subroutine atmos_global_diag_init (axes, area) +integer, intent(in), dimension(:) :: axes +real, intent(in), dimension(:,:) :: area + +!----------------------------------------------------------------------- +! local data + +integer :: iunit, ierr, io + +!----------------------------------------------------------------------- + + if (module_is_initialized) then + call error_mesg ('atmos_global_diag_mod', & + 'module has already been initialized', WARNING) + return + endif + +!----------------------------------------------------------------------- +!----- read namelist ----- + read (input_nml_file, nml=atmos_global_diag_nml, iostat=io) + ierr = check_nml_error (io, 'atmos_global_diag_nml') + +!----- write version and namelist to log file ----- + + iunit = stdlog() + call write_version_number ( version, tagname ) + if (mpp_pe() == mpp_root_pe()) write (iunit, nml=atmos_global_diag_nml) + +!----------------------------------------------------------------------- +!----- allocate storage for fields ----- + allocate(fields(max_fields)) + num_fields = 0 + +!----- determine the atmospheric 2D domain from the axes using diag_manager ----- + Domain2 = get_domain2d(axes) + if (Domain2 .eq. null_domain2d) then + call error_mesg ('atmos_global_diag_mod', & + 'could not determine 2d domain', FATAL) + endif + +!----- save cell areas and the computed global area ---- + id = size(area,1); jd = size(area,2) + allocate(area_g(id,jd)) + area_g = area + area_g_sum = mpp_global_sum (Domain2, area_g, flags=BITWISE_EFP_SUM) + +! null axis for global fields in diag_manager +! null_axis(1) = null_axis_id + + module_is_initialized = .true. + +!----------------------------------------------------------------------- + +end subroutine atmos_global_diag_init + +!####################################################################### + +function register_global_diag_field (field_name, Time_init, long_name, & + units, standard_name, realm, buffer, use_masking ) +character(len=*), intent(in) :: field_name +type(time_type), intent(in) :: Time_init +character(len=*), intent(in), optional :: long_name, units, standard_name, realm +logical, intent(in), optional :: buffer, use_masking + +!----------------------------------------------------------------------- + +integer :: register_global_diag_field +integer :: field_index, field_id + +!----------------------------------------------------------------------- + + if (.not.module_is_initialized) call error_mesg & + ('atmos_global_diag_mod', 'module has not been initialized', FATAL) + + register_global_diag_field = 0 + + field_index = find_field_index(field_name) + + if (field_index > 0) call error_mesg ('atmos_global_diag_mod', & + 'trying to register field already registered', FATAL) + + ! register field with diag_manager + field_id = register_diag_field (mod_name, field_name, Time_init, & + long_name, units, standard_name=standard_name, & + realm=realm) + if (field_id == DIAG_FIELD_NOT_FOUND) return + + ! register field with this module + num_fields = num_fields+1 + if (num_fields > max_fields) call error_mesg ('atmos_global_diag_mod', & + 'exceeded max_fields, increase value in namelist', FATAL) + register_global_diag_field = num_fields + + fields(num_fields)%field_name = trim(field_name) + fields(num_fields)%field_id = field_id + fields(num_fields)%use_masking = .false. + fields(num_fields)%use_buffer = .false. + + if (present(buffer)) then + if (buffer) then + allocate(fields(num_fields)%buffer(id,jd)) + fields(num_fields)%buffer = missing_value + fields(num_fields)%use_buffer = .true. + if (present(use_masking)) fields(num_fields)%use_masking = use_masking + endif + endif + + return + +!----------------------------------------------------------------------- + +end function register_global_diag_field + +!####################################################################### + +function get_global_diag_field_id (field_num) +integer, intent(in) :: field_num + +integer :: get_global_diag_field_id +!----------------------------------------------------------------------- + get_global_diag_field_id = DIAG_FIELD_NOT_FOUND + + if (.not.module_is_initialized) call error_mesg & + ('atmos_global_diag_mod', 'module has not been initialized', FATAL) + + if (field_num == 0) return + if (field_num < 0 .or. field_num > num_fields) call error_mesg & + ('atmos_global_diag_mod', 'invalid field number in get_global_diag_field_id', FATAL) + + get_global_diag_field_id = fields(field_num)%field_id + +!----------------------------------------------------------------------- + +end function get_global_diag_field_id + +!####################################################################### + +subroutine buffer_global_diag (field_num, data, Time, is_in, js_in, mask) +integer, intent(in) :: field_num +real, intent(in) :: data(:,:) +type(time_type), intent(in) :: Time +integer, optional, intent(in) :: is_in, js_in +logical, optional, intent(in) :: mask(:,:) +!----------------------------------------------------------------------- + +integer :: is, ie, js, je + +!----------------------------------------------------------------------- + + if (.not.module_is_initialized) call error_mesg & + ('atmos_global_diag_mod', 'module has not been initialized', FATAL) + + if (field_num == 0) return + if (field_num < 0 .or. field_num > num_fields) call error_mesg & + ('atmos_global_diag_mod', 'invalid field number in buffer_global_data', FATAL) + + if (.not.fields(field_num)%use_buffer) call error_mesg ('atmos_global_diag_mod', & + 'buffer not allocated in buffer_global_data for field='//trim(fields(field_num)%field_name), FATAL) + + is = 1; if (present(is_in)) is = is_in + js = 1; if (present(js_in)) js = js_in + ie = is + size(data,1) -1 + je = js + size(data,2) -1 + + if (present(mask)) then + if (.not.fields(field_num)%use_masking) call error_mesg & + ('atmos_global_diag_mod', 'attempt to pass mask field without '// & + 'setting use_masking=True', FATAL) + where (mask) fields(field_num)%buffer(is:ie,js:je) = data + else + fields(field_num)%buffer(is:ie,js:je) = data + endif + fields(field_num)%Time = Time + +!----------------------------------------------------------------------- + +end subroutine buffer_global_diag + +!####################################################################### + +function send_global_diag_data (field_num, data, Time, area, mask) +integer, intent(in) :: field_num +real, intent(in) :: data(:,:) +type(time_type), intent(in) :: Time +real, optional, intent(in) :: area(:,:) +logical, optional, intent(in) :: mask(:,:) + +logical :: send_global_diag_data +!----------------------------------------------------------------------- + +real, dimension(id,jd) :: data_gbl, area_gbl +real :: gbl_sum, area_sum + +!----------------------------------------------------------------------- + + if (.not.module_is_initialized) call error_mesg & + ('atmos_global_diag_mod', 'module has not been initialized', FATAL) + + if (field_num == 0) return + + if (field_num < 0 .or. field_num > num_fields) call error_mesg & + ('atmos_global_diag_mod', 'invalid field number in send_global_diag_data', FATAL) + + if (fields(field_num)%use_buffer) call error_mesg ('atmos_global_diag_mod', & + 'buffer allocated in send_global_diag_data for field='//trim(fields(field_num)%field_name), FATAL) + + if (.not.present(area) .and. .not.present(mask)) then + data_gbl = data*area_g + area_sum = area_g_sum + + else if (present(area) .and. present(mask)) then + where (mask) + data_gbl = data*area + area_gbl = area + elsewhere + data_gbl = 0.0 + area_gbl = 0.0 + endwhere + area_sum = mpp_global_sum (Domain2, area_gbl, flags=BITWISE_EFP_SUM) + + else if (present(area) .and. .not.present(mask)) then + data_gbl = data*area + area_sum = mpp_global_sum (Domain2, area, flags=BITWISE_EFP_SUM) + + else if (.not.present(area) .and. present(mask)) then + where (mask) + data_gbl = data + area_gbl = area_g + elsewhere + data_gbl = 0.0 + area_gbl = 0.0 + endwhere + area_sum = mpp_global_sum (Domain2, area_gbl, flags=BITWISE_EFP_SUM) + endif + + gbl_sum = mpp_global_sum (Domain2, data_gbl, flags=BITWISE_EFP_SUM) + send_global_diag_data = send_data (fields(field_num)%field_id, gbl_sum/area_sum, Time) + +!----------------------------------------------------------------------- + +end function send_global_diag_data + +!####################################################################### + +function send_global_diag_buffer (field_num) +integer, intent(in) :: field_num +logical :: send_global_diag_buffer +!----------------------------------------------------------------------- + +real, dimension(id,jd) :: data_g, area_m +real :: gbl_sum, area_sum + +!----------------------------------------------------------------------- + + if (.not.module_is_initialized) call error_mesg & + ('atmos_global_diag_mod', 'module has not been initialized', FATAL) + + if (field_num == 0) return + + if (field_num < 0 .or. field_num > num_fields) call error_mesg & + ('atmos_global_diag_mod', 'invalid field number in send_global_diag_buffer', FATAL) + + if (.not.fields(field_num)%use_buffer) call error_mesg ('atmos_global_diag_mod', & + 'buffer not allocated in send_global_diag_buffer for field='//trim(fields(field_num)%field_name), FATAL) + + ! weight data with area + data_g = fields(field_num)%buffer*area_g + area_sum = area_g_sum + + !--- when masking cells, need to compute new global area --- + if (fields(field_num)%use_masking) then + where(fields(field_num)%buffer .eq. missing_value) + data_g = 0.0 + area_m = 0.0 + elsewhere + area_m = area_g + endwhere + ! global sum of masked area + area_sum = mpp_global_sum(Domain2, area_m, flags=BITWISE_EFP_SUM) + endif + + ! sum of data*area + gbl_sum = mpp_global_sum (Domain2, data_g, flags=BITWISE_EFP_SUM) + + send_global_diag_buffer = send_data (fields(field_num)%field_id, gbl_sum/area_sum, & + fields(field_num)%Time) + +!----------------------------------------------------------------------- + +end function send_global_diag_buffer + +!####################################################################### + +subroutine atmos_global_diag_end + +integer :: ind + +!----------------------------------------------------------------------- + + if (.not.module_is_initialized) then + call error_mesg ('atmos_global_diag_mod', & + 'module has not been initialized, '// & + 'nothing to do in atmos_global_diag_end', WARNING) + return + endif + + do ind = 1, num_fields + if (fields(ind)%use_buffer) then + deallocate(fields(ind)%buffer) + fields(ind)%use_buffer = .false. + endif + enddo + deallocate(fields) + + module_is_initialized = .false. + +!----------------------------------------------------------------------- + +end subroutine atmos_global_diag_end + +!####################################################################### +! private routine + +function find_field_index (field_name) +character(len=*), intent(in) :: field_name +integer :: find_field_index +integer :: ind + + find_field_index = 0 + + do ind = 1, num_fields + if (trim(field_name) .eq. trim(fields(ind)%field_name)) then + find_field_index = ind + return + endif + enddo + +end function find_field_index + +!####################################################################### + +end module atmos_global_diag_mod + diff --git a/atmos_shared/atmos_tracer_driver.F90 b/atmos_shared/atmos_tracer_driver.F90 new file mode 100644 index 00000000..2ad9fb1e --- /dev/null +++ b/atmos_shared/atmos_tracer_driver.F90 @@ -0,0 +1,427 @@ +!*********************************************************************** +!* GNU Lesser General Public License +!* +!* This file is part of the GFDL Atmosphere Null Model Component. +!* +!* Atmos Null is free software: you can redistribute it and/or modify it +!* under the terms of the GNU Lesser General Public License as published +!* by the Free Software Foundation, either version 3 of the License, or +!* (at your option) any later version. +!* +!* Atmos Null is distributed in the hope that it will be useful, but +!* WITHOUT ANY WARRANTY; without even the implied warranty of +!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +!* General Public License for more details. +!* +!* You should have received a copy of the GNU Lesser General Public +!* License along with Atmos Null. +!* If not, see . +!*********************************************************************** + +module atmos_tracer_driver_mod +! +! William Cooke +! + +! +! Matt Harrison +! + +! +! Bruce Wyman +! + +! + +! +! This code allows the user to easily add tracers to the FMS framework. +! + +! +! +! This is a null module to be used when no atmospheric model is required. +! +! This code allows a user to easily implement tracer code in the FMS +! framework. The tracer and tracer tendency arrays are supplied along +! with longtitude, latitude, wind, temperature, and pressure +! information which allows a user to implement sources and sinks of the +! tracer which depend on these parameters. +! +! In the following example, radon being implemented in the atmosphere +! will be used as an example of how to implement a tracer in the FMS +! framework. +! +! Within the global scope of tracer_driver_mod a use +! statement should be inserted for each tracer to be added. +!
      use radon_mod, only : radon_sourcesink, radon_init, radon_end
+! +! An integer parameter, which will be used as an identifier for the +! tracer, should be assigned. +!
+!      integer :: nradon 
+!
+! Within tracer_driver_init a call to the tracer manager is needed in +! order to identify which tracer it has set the tracer as. +!
+!      nradon = get_tracer_index(MODEL_ATMOS,'radon')
+!
+! Here MODEL_ATMOS is a parameter defined in field_manager. +! 'radon' is the name of the tracer within the field_table. +! +! If the tracer exists then the integer returned will be positive and it +! can be used to call the initialization routines for the individual +! tracers. +!
+!      if (nradon > 0) then
+!           call radon_init(Argument list)
+!      endif
+!
+! +! Within tracer_driver the user can also use the identifier to surround +! calls to the source-sink routines for the tracer of interest. +! +!
+!      if (nradon > 0 .and. nradon <= nt) then
+!          call radon_sourcesink (Argument list)
+!          rdt(:,:,:,nradon)=rdt(:,:,:,nradon)+rtnd(:,:,:)
+!      endif
+!
+! +! It is the users responsibility to add the tendency generated by the +! sourcesink routine. +! +! Within tracer_driver_end the user can add calls to the +! terminators for the appropriate source sink routines. +! +!
      call radon_end
+! +! This may simply be a deallocation statement or a routine to send +! output to the logfile stating that the termination routine has been +! called. +! + + +!----------------------------------------------------------------------- + +use time_manager_mod, only : time_type + +implicit none +private +!----------------------------------------------------------------------- +!----- interfaces ------- + +public atmos_tracer_driver, atmos_tracer_flux_init +public atmos_tracer_driver_init, atmos_tracer_driver_end +public atmos_tracer_driver_gather_data +public atmos_tracer_driver_gather_data_down +public atmos_tracer_has_surf_setl_flux, get_atmos_tracer_surf_setl_flux + +!----------------------------------------------------------------------- +!----------- namelist ------------------- +!----------------------------------------------------------------------- +! +! When initializing additional tracers, the user needs to make the +! following changes. +! +! Add an integer variable below for each additional tracer. +! This should be initialized to zero. +! +!----------------------------------------------------------------------- + +character(len=6), parameter :: module_name = 'tracer' + +logical :: module_is_initialized = .FALSE. + + +!---- version number ----- +character(len=128) :: version = '$Id$' +character(len=128) :: tagname = '$Name$' +!----------------------------------------------------------------------- + +contains + +!####################################################################### + +! +! +! A routine which allows tracer code to be called. +! +! +! This subroutine calls the source sink routines for atmospheric +! tracers. This is the interface between the dynamical core of the +! model and the tracer code. It should supply all the necessary +! information to a user that they need in order to calculate the +! tendency of that tracer with respect to emissions or chemical losses. +! +! +! +! +! Local domain boundaries. +! +! +! Model time. +! +! +! Longitude of the centre of the model gridcells +! +! +! Latitude of the centre of the model gridcells +! +! +! Land/sea mask. +! +! +! Pressures on the model half levels. +! +! +! Pressures on the model full levels. +! +! +! The tracer array in the component model. +! +! +! Zonal wind speed. +! +! +! Meridonal wind speed. +! +! +! Temperature. +! +! +! Specific humidity. This may also be accessible as a +! portion of the tracer array. +! +! +! Friction velocity :: +! The magnitude of the wind stress is density*(ustar**2) +! The drag coefficient for momentum is u_star**2/(u**2+v**2) +! +! +! The tendency of the tracer array in the compenent +! model. The tendency due to sources and sinks computed +! in the individual tracer routines should be added to +! this array before exiting tracer_driver. +! +! +! The tracer array in the component model for the previous timestep. +! +! +! The array of diagnostic tracers. As these may be changed within the +! tracer routines for diagnostic purposes, they need to be writable. +! +! +! Integer array describing which model layer intercepts the surface. +! + subroutine atmos_tracer_driver (is, ie, js, je, Time, lon, lat, & + area, z_pbl, rough_mom, & + land, phalf, pfull, & + u, v, t, q, r, & + rm, rdt, dt, & + u_star, b_star, q_star, & + z_half, z_full,& + t_surf_rad, albedo, & + Time_next, mask, & + kbot) + +!----------------------------------------------------------------------- +integer, intent(in) :: is, ie, js, je +type(time_type), intent(in) :: Time +real, intent(in), dimension(:,:) :: lon, lat +real, intent(in), dimension(:,:) :: u_star, b_star, q_star +real, intent(in), dimension(:,:) :: land +real, intent(in), dimension(:,:) :: area, z_pbl, rough_mom +real, intent(in), dimension(:,:,:) :: phalf, pfull +real, intent(in), dimension(:,:,:) :: u, v, t, q +real, intent(inout), dimension(:,:,:,:) :: r +real, intent(inout), dimension(:,:,:,:) :: rm +real, intent(inout), dimension(:,:,:,:) :: rdt +real, intent(in) :: dt !timestep(used in chem_interface) +real, intent(in), dimension(:,:,:) :: z_half !height in meters at half levels +real, intent(in), dimension(:,:,:) :: z_full !height in meters at full levels +real, intent(in), dimension(:,:) :: t_surf_rad !surface temperature +real, intent(in), dimension(:,:) :: albedo +type(time_type), intent(in) :: Time_next +integer, intent(in), dimension(:,:), optional :: kbot +real, intent(in), dimension(:,:,:), optional :: mask + + + return + + end subroutine atmos_tracer_driver +! + +!####################################################################### + +! +! +! Subroutine to initialize the ocean-atmosphere gas flux modules +! +! +! Subroutine to initialize the ocean-atmosphere gas flux modules +! + +subroutine atmos_tracer_flux_init + +return + +end subroutine atmos_tracer_flux_init +! + +!####################################################################### + +! +! +! Subroutine to initialize the tracer driver module. +! +! +! The purpose of the arguments here are for passing on to the individual +! tracer code. The user may wish to provide initial values which can be +! implemented in the initialization part of the tracer code. Remember that +! the tracer manager will provide a simple fixed or exponential profile if +! the user provides data for this within the field table. However if a more +! complicated profile is required then it should be set up in the +! initialization section of the user tracer code. +! +! +! +! The longitudes for the local domain. +! +! +! The latitudes for the local domain. +! +! +! optional mask (0. or 1.) that designates which grid points +! are above (=1.) or below (=0.) the ground dimensioned as +! (nlon,nlat,nlev). +! +! +! Model time. +! +! +! The axes relating to the tracer array dimensioned as +! (nlon, nlat, nlev, ntime) +! +! +! Tracer fields dimensioned as (nlon,nlat,nlev,ntrace). +! + subroutine atmos_tracer_driver_init (lonb, latb, r, axes, Time, phalf, mask) + +!----------------------------------------------------------------------- + real, intent(in), dimension(:,:) :: lonb, latb + real, intent(inout), dimension(:,:,:,:) :: r +type(time_type), intent(in) :: Time + integer, intent(in) :: axes(4) + real, intent(in), dimension(:,:,:) :: phalf + real, intent(in), dimension(:,:,:), optional :: mask + +!----------------------------------------------------------------------- +! +! When initializing additional tracers, the user needs to make changes +! +!----------------------------------------------------------------------- + + module_is_initialized = .TRUE. + + return + + end subroutine atmos_tracer_driver_init +! + +!####################################################################### + +! +! +! Subroutine to terminate the tracer driver module. +! +! +! Termination routine for tracer_driver. It should also call +! the destructors for the individual tracer routines. +! +! + subroutine atmos_tracer_driver_end + +!----------------------------------------------------------------------- + + module_is_initialized = .FALSE. + +!----------------------------------------------------------------------- + + return + + end subroutine atmos_tracer_driver_end +! + +!####################################################################### + +! +! +! Subroutine to terminate the tracer driver module. +! +! +! Termination routine for tracer_driver. It should also call +! the destructors for the individual tracer routines. +! +! + subroutine atmos_tracer_driver_gather_data(gas_fields, tr_bot) + +use coupler_types_mod, only: coupler_2d_bc_type + +type(coupler_2d_bc_type), intent(inout) :: gas_fields +real, dimension(:,:,:), intent(in) :: tr_bot + +!----------------------------------------------------------------------- + +!----------------------------------------------------------------------- + + return + + end subroutine atmos_tracer_driver_gather_data + +!###################################################################### + subroutine atmos_tracer_driver_gather_data_down(gas_fields, tr_bot) + + use coupler_types_mod, only: coupler_2d_bc_type + + type(coupler_2d_bc_type), intent(inout) :: gas_fields + real, dimension(:,:,:), intent(in) :: tr_bot + + return + + end subroutine atmos_tracer_driver_gather_data_down +! +!###################################################################### +! given a tracer index, returns true if this tracer has non-zero +! sedimentation flux at the bottom of the atmosphere +function atmos_tracer_has_surf_setl_flux(tr) result(ret) + logical :: ret + integer, intent(in) :: tr ! tracer index + + ret=.FALSE. +end function + +!###################################################################### +! given a tracer index, returns sedimentation flux at the bottom of +! the atmosphere for this tracer +subroutine get_atmos_tracer_surf_setl_flux(tr, setl_flux, dsetl_dtr) + integer, intent(in) :: tr ! tracer index + real, intent(out) :: setl_flux(:,:) ! sedimentation flux at the bottom of the atmosphere + real, intent(out) :: dsetl_dtr(:,:) ! derivative of sedimentation flux w.r.t. + ! the tracer concentration in the bottom layer + + setl_flux(:,:) = 0.0 ; dsetl_dtr(:,:) = 0.0 + return +end subroutine +!###################################################################### + + +end module atmos_tracer_driver_mod