!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS 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 General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module ice_bergs use constants_mod, only: radius, pi, omega, HLF use fms_mod, only: open_namelist_file, check_nml_error, close_file use fms_mod, only: field_exist, get_global_att_value use fms_mod, only: stdlog, stderr, error_mesg, FATAL, WARNING use fms_mod, only: write_version_number, read_data, write_data, file_exist use mosaic_mod, only: get_mosaic_ntiles, get_mosaic_ncontacts use mpp_mod, only: mpp_pe, mpp_root_pe, mpp_sum, mpp_min, mpp_max, NULL_PE use mpp_mod, only: mpp_send, mpp_recv, mpp_sync_self, mpp_chksum use mpp_mod, only: mpp_clock_begin, mpp_clock_end, mpp_clock_id use mpp_mod, only: CLOCK_COMPONENT, CLOCK_SUBCOMPONENT, CLOCK_LOOP use fms_mod, only: clock_flag_default use mpp_domains_mod, only: domain2D, mpp_update_domains, mpp_define_domains use mpp_domains_mod, only: SCALAR_PAIR, CGRID_NE, BGRID_NE, CORNER use mpp_domains_mod, only: mpp_get_compute_domain, mpp_get_data_domain use mpp_domains_mod, only: CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE use mpp_domains_mod, only: mpp_get_neighbor_pe, NORTH, SOUTH, EAST, WEST use mpp_domains_mod, only: mpp_define_io_domain use time_manager_mod, only: time_type, get_date, get_time, set_date, operator(-) use diag_manager_mod, only: register_diag_field, register_static_field, send_data use diag_manager_mod, only: diag_axis_init implicit none ; private include 'netcdf.inc' public icebergs_init, icebergs_end, icebergs_run, icebergs_stock_pe public icebergs_incr_mass, icebergs_save_restart type :: icebergs_gridded type(domain2D), pointer :: domain ! MPP domain integer :: halo ! Nominal halo width integer :: isc, iec, jsc, jec ! Indices of computational domain integer :: isd, ied, jsd, jed ! Indices of data domain integer :: my_pe, pe_N, pe_S, pe_E, pe_W ! MPI PE identifiers real, dimension(:,:), pointer :: lon=>NULL() ! Longitude of cell corners real, dimension(:,:), pointer :: lat=>NULL() ! Latitude of cell corners real, dimension(:,:), pointer :: lonc=>NULL() ! Longitude of cell centers real, dimension(:,:), pointer :: latc=>NULL() ! Latitude of cell centers real, dimension(:,:), pointer :: dx=>NULL() ! Length of cell edge (m) real, dimension(:,:), pointer :: dy=>NULL() ! Length of cell edge (m) real, dimension(:,:), pointer :: area=>NULL() ! Area of cell (m^2) real, dimension(:,:), pointer :: msk=>NULL() ! Ocean-land mask (1=ocean) real, dimension(:,:), pointer :: cos=>NULL() ! Cosine from rotation matrix to lat-lon coords real, dimension(:,:), pointer :: sin=>NULL() ! Sine from rotation matrix to lat-lon coords real, dimension(:,:), pointer :: uo=>NULL() ! Ocean zonal flow (m/s) real, dimension(:,:), pointer :: vo=>NULL() ! Ocean meridional flow (m/s) real, dimension(:,:), pointer :: ui=>NULL() ! Ice zonal flow (m/s) real, dimension(:,:), pointer :: vi=>NULL() ! Ice meridional flow (m/s) real, dimension(:,:), pointer :: ua=>NULL() ! Atmosphere zonal flow (m/s) real, dimension(:,:), pointer :: va=>NULL() ! Atmosphere meridional flow (m/s) real, dimension(:,:), pointer :: ssh=>NULL() ! Sea surface height (m) real, dimension(:,:), pointer :: sst=>NULL() ! Sea surface temperature (oC) real, dimension(:,:), pointer :: cn=>NULL() ! Sea-ice concentration (0 to 1) real, dimension(:,:), pointer :: hi=>NULL() ! Sea-ice thickness (m) real, dimension(:,:), pointer :: calving=>NULL() ! Calving mass rate [frozen runoff] (kg/s) (into stored ice) real, dimension(:,:), pointer :: calving_hflx=>NULL() ! Calving heat flux [heat content of calving] (W/m2) (into stored ice) real, dimension(:,:), pointer :: floating_melt=>NULL() ! Net melting rate to icebergs + bits (kg/s/m^2) real, dimension(:,:), pointer :: berg_melt=>NULL() ! Melting+erosion rate of icebergs (kg/s/m^2) real, dimension(:,:), pointer :: melt_buoy=>NULL() ! Buoyancy componenet of melting rate (kg/s/m^2) real, dimension(:,:), pointer :: melt_eros=>NULL() ! Erosion component of melting rate (kg/s/m^2) real, dimension(:,:), pointer :: melt_conv=>NULL() ! Convective component of melting rate (kg/s/m^2) real, dimension(:,:), pointer :: bergy_src=>NULL() ! Mass flux from berg erosion into bergy bits (kg/s/m^2) real, dimension(:,:), pointer :: bergy_melt=>NULL() ! Melting rate of bergy bits (kg/s/m^2) real, dimension(:,:), pointer :: bergy_mass=>NULL() ! Mass distribution of bergy bits (kg/s/m^2) real, dimension(:,:), pointer :: virtual_area=>NULL() ! Virtual surface coverage by icebergs (m^2) real, dimension(:,:), pointer :: mass=>NULL() ! Mass distribution (kg/m^2) real, dimension(:,:,:), pointer :: mass_on_ocean=>NULL() ! Mass distribution partitioned by neighbor (kg/m^2) real, dimension(:,:), pointer :: tmp=>NULL() ! Temporary work space real, dimension(:,:), pointer :: tmpc=>NULL() ! Temporary work space real, dimension(:,:,:), pointer :: stored_ice=>NULL() ! Accumulated ice mass flux at calving locations (kg) real, dimension(:,:), pointer :: stored_heat=>NULL() ! Heat content of stored ice (J) real, dimension(:,:,:), pointer :: real_calving=>NULL() ! Calving rate into iceberg class at calving locations (kg/s) real, dimension(:,:), pointer :: iceberg_heat_content=>NULL() ! Distributed heat content of bergs (J/m^2) ! Diagnostics handles integer :: id_uo=-1, id_vo=-1, id_calving=-1, id_stored_ice=-1, id_accum=-1, id_unused=-1, id_floating_melt=-1 integer :: id_melt_buoy=-1, id_melt_eros=-1, id_melt_conv=-1, id_virtual_area=-1, id_real_calving=-1 integer :: id_calving_hflx_in=-1, id_stored_heat=-1, id_melt_hflx=-1, id_heat_content=-1 integer :: id_mass=-1, id_ui=-1, id_vi=-1, id_ua=-1, id_va=-1, id_sst=-1, id_cn=-1, id_hi=-1 integer :: id_bergy_src=-1, id_bergy_melt=-1, id_bergy_mass=-1, id_berg_melt=-1 end type icebergs_gridded type :: xyt real :: lon, lat, day real :: mass, thickness, width, length, uvel, vvel real :: uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi real :: mass_of_bits, heat_density integer :: year type(xyt), pointer :: next=>NULL() end type xyt type :: iceberg type(iceberg), pointer :: prev=>NULL(), next=>NULL() ! State variables (specific to the iceberg, needed for restarts) real :: lon, lat, uvel, vvel, mass, thickness, width, length real :: start_lon, start_lat, start_day, start_mass, mass_scaling real :: mass_of_bits, heat_density integer :: start_year integer :: ine, jne ! nearest index in NE direction (for convenience) real :: xi, yj ! Non-dimensional coords within current cell (0..1) ! Environment variables (as seen by the iceberg) real :: uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi type(xyt), pointer :: trajectory=>NULL() end type iceberg type :: buffer integer :: size=0 real, dimension(:,:), pointer :: data end type buffer type, public :: icebergs ; private type(icebergs_gridded), pointer :: grd type(iceberg), pointer :: first=>NULL() type(xyt), pointer :: trajectories=>NULL() real :: dt ! Time-step between iceberg calls (should make adaptive?) integer :: current_year real :: current_yearday ! 1.00-365.99 integer :: traj_sample_hrs integer :: verbose_hrs integer :: clock, clock_mom, clock_the, clock_int, clock_cal, clock_com, clock_ini, clock_ior, clock_iow, clock_dia ! ids for fms timers real :: rho_bergs ! Density of icebergs real :: LoW_ratio ! Initial ratio L/W for newly calved icebergs real :: bergy_bit_erosion_fraction ! Fraction of erosion melt flux to divert to bergy bits real :: sicn_shift ! Shift of sea-ice concentration in erosion flux modulation (0NULL(), ibuffer_n=>NULL() type(buffer), pointer :: obuffer_s=>NULL(), ibuffer_s=>NULL() type(buffer), pointer :: obuffer_e=>NULL(), ibuffer_e=>NULL() type(buffer), pointer :: obuffer_w=>NULL(), ibuffer_w=>NULL() ! Budgets real :: net_calving_received=0., net_calving_returned=0. real :: net_incoming_calving=0., net_outgoing_calving=0. real :: net_incoming_calving_heat=0., net_outgoing_calving_heat=0. real :: net_incoming_calving_heat_used=0., net_heat_to_bergs=0. real :: stored_start=0., stored_end=0. real :: stored_heat_start=0., stored_heat_end=0., net_heat_to_ocean=0. real :: net_calving_used=0., net_calving_to_bergs=0. real :: floating_mass_start=0., floating_mass_end=0. real :: floating_heat_start=0., floating_heat_end=0. real :: icebergs_mass_start=0., icebergs_mass_end=0. real :: bergy_mass_start=0., bergy_mass_end=0. real :: returned_mass_on_ocean=0. real :: net_melt=0., berg_melt=0., bergy_src=0., bergy_melt=0. integer :: nbergs_calved=0, nbergs_melted=0, nbergs_start=0, nbergs_end=0 integer :: nspeeding_tickets=0 integer, dimension(:), pointer :: nbergs_calved_by_class=>NULL() end type icebergs ! Global constants character(len=*), parameter :: version = '$Id: ice_bergs.F90,v 1.1.2.113.2.1 2009/10/23 17:48:33 nnz Exp $' character(len=*), parameter :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' integer, parameter :: nclasses=10 ! Number of ice bergs classes integer, parameter :: file_format_major_version=0 integer, parameter :: file_format_minor_version=1 integer, parameter :: delta_buf=25 ! Size by which to increment buffers real, parameter :: pi_180=pi/180. ! Converts degrees to radians real, parameter :: rho_ice=916.7 ! Density of fresh ice @ 0oC (kg/m^3) real, parameter :: rho_water=999.8 ! Density of fresh water @ 0oC (kg/m^3) real, parameter :: rho_air=1.1 ! Density of air @ 0oC (kg/m^3) ??? real, parameter :: rho_seawater=1025. ! Approx. density of surface sea water @ 0oC (kg/m^3) real, parameter :: gravity=9.8 ! Gravitational acceleratio (m/s^2) real, parameter :: Cd_av=1.3 ! (Vertical) Drag coefficient between bergs and atmos (?) real, parameter :: Cd_ah=0.0055 ! (Horizontal) Drag coefficient between bergs and atmos (?) real, parameter :: Cd_wv=0.9 ! (Vertical) Drag coefficient between bergs and ocean (?) real, parameter :: Cd_wh=0.0012 ! (Horizontal) Drag coefficient between bergs and ocean (?) real, parameter :: Cd_iv=0.9 ! (Vertical) Drag coefficient between bergs and sea-ice (?) !TOM> no horizontal drag for sea ice! real, parameter :: Cd_ih=0.0012 ! (Horizontal) Drag coefficient between bergs and sea-ice (?) ! Global data (minimal for debugging) logical :: verbose=.false. ! Be verbose to stderr logical :: budget=.true. ! Calculate budgets logical :: debug=.false. ! Turn on debugging logical :: really_debug=.false. ! Turn on debugging logical :: parallel_reprod=.true. ! Reproduce across different PE decompositions logical :: use_slow_find=.true. ! Use really slow (but robust) find_cell for reading restarts logical :: ignore_ij_restart=.false. ! Read i,j location from restart if available (needed to use restarts on different grids) logical :: generate_test_icebergs=.false. ! Create icebergs in absence of a restart file contains ! ############################################################################## subroutine accel(bergs, berg, i, j, xi, yj, lat, uvel, vvel, uvel0, vvel0, dt, ax, ay, debug_flag) ! Arguments type(icebergs), pointer :: bergs type(iceberg), pointer :: berg integer, intent(in) :: i, j real, intent(in) :: xi, yj, lat, uvel, vvel, uvel0, vvel0, dt real, intent(inout) :: ax, ay logical, optional :: debug_flag ! Local variables type(icebergs_gridded), pointer :: grd real :: uo, vo, ui, vi, ua, va, uwave, vwave, ssh_x, ssh_y, sst, cn, hi real :: f_cori, T, D, W, L, M, F real :: drag_ocn, drag_atm, drag_ice, wave_rad real :: c_ocn, c_atm, c_ice real :: ampl, wmod, Cr, Lwavelength, Lcutoff, Ltop real, parameter :: alpha=0.0, beta=1.0, accel_lim=1.e-2, Cr0=0.06, vel_lim=15. real :: lambda, detA, A11, A12, axe, aye, D_hi real :: uveln, vveln, us, vs, speed, loc_dx, new_speed logical :: dumpit integer :: itloop ! For convenience grd=>bergs%grd ! Interpolate gridded fields to berg call interp_flds(grd, i, j, xi, yj, uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi) f_cori=(2.*omega)*sin(pi_180*lat) M=berg%mass T=berg%thickness ! total thickness D=(bergs%rho_bergs/rho_seawater)*T ! draught (keel depth) F=T-D ! freeboard W=berg%width L=berg%length hi=min(hi,D) D_hi=max(0.,D-hi) ! Wave radiation uwave=ua-uo; vwave=va-vo ! Use wind speed rel. to ocean for wave model (aja)? wmod=uwave*uwave+vwave*vwave ! The wave amplitude and length depend on the wind speed relative to the ocean current; ! actually wmod is wmod**2 here. ampl=0.5*0.02025*wmod ! This is "a", the wave amplitude Lwavelength=0.32*wmod ! Surface wave length fitted to data in table at ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html Lcutoff=0.125*Lwavelength Ltop=0.25*Lwavelength Cr=Cr0*min(max(0.,(L-Lcutoff)/((Ltop-Lcutoff)+1.e-30)),1.) ! Wave radiation coefficient ! fitted to graph from Carrieres et al., POAC Drift Model. wave_rad=0.5*rho_seawater/M*Cr*gravity*ampl*min(ampl,F)*(2.*W*L)/(W+L) wmod = sqrt(ua*ua+va*va) ! Wind speed if (wmod.ne.0.) then uwave=ua/wmod ! Wave radiation force acts in wind direction ... vwave=va/wmod else uwave=0.; vwave=0.; wave_rad=0. ! ... and only when wind is present. endif ! Weighted drag coefficients c_ocn=rho_seawater/M*(0.5*Cd_wv*W*(D_hi)+Cd_wh*W*L) c_atm=rho_air /M*(0.5*Cd_av*W*F +Cd_ah*W*L) c_ice=rho_ice /M*(0.5*Cd_iv*W*hi ) if (abs(ui)+abs(vi).eq.0.) c_ice=0. uveln=uvel; vveln=vvel ! Copy starting uvel, vvel do itloop=1,2 ! Iterate on drag coefficients us=0.5*(uveln+uvel); vs=0.5*(vveln+vvel) drag_ocn=c_ocn*sqrt( (us-uo)**2+(vs-vo)**2 ) drag_atm=c_atm*sqrt( (us-ua)**2+(vs-va)**2 ) drag_ice=c_ice*sqrt( (us-ui)**2+(vs-vi)**2 ) ! Explicit accelerations !axe= f_cori*vvel -gravity*ssh_x +wave_rad*uwave & ! -drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui) !aye=-f_cori*uvel -gravity*ssh_y +wave_rad*vwave & ! -drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi) axe=-gravity*ssh_x +wave_rad*uwave aye=-gravity*ssh_y +wave_rad*vwave if (alpha>0.) then ! If implicit, use time-level (n) rather than RK4 latest axe=axe+f_cori*vvel0 aye=aye-f_cori*uvel0 else axe=axe+f_cori*vvel aye=aye-f_cori*uvel endif if (beta>0.) then ! If implicit, use time-level (n) rather than RK4 latest axe=axe-drag_ocn*(uvel0-uo) -drag_atm*(uvel0-ua) -drag_ice*(uvel0-ui) aye=aye-drag_ocn*(vvel0-vo) -drag_atm*(vvel0-va) -drag_ice*(vvel0-vi) else axe=axe-drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui) aye=aye-drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi) endif ! Solve for implicit accelerations if (alpha+beta.gt.0.) then lambda=drag_ocn+drag_atm+drag_ice A11=1.+beta*dt*lambda A12=alpha*dt*f_cori detA=1./(A11**2+A12**2) ax=detA*(A11*axe+A12*aye) ay=detA*(A11*aye-A12*axe) else ax=axe; ay=aye endif uveln=uvel0+dt*ax vveln=vvel0+dt*ay enddo ! itloop ! Limit speed of bergs based on a CFL criteria if (bergs%speed_limit>0.) then speed=sqrt(uveln*uveln+vveln*vveln) ! Speed of berg if (speed>0.) then loc_dx=min(0.5*(grd%dx(i,j)+grd%dx(i,j-1)),0.5*(grd%dy(i,j)+grd%dy(i-1,j))) ! min(dx,dy) !new_speed=min(loc_dx/dt*bergs%speed_limit,speed) ! Restrict speed to dx/dt x factor new_speed=loc_dx/dt*bergs%speed_limit ! Speed limit as a factor of dx / dt if (new_speedvel_lim.or.abs(vveln)>vel_lim) then dumpit=.true. write(stderr(),'("pe=",i3,x,a)') mpp_pe(),'Dump triggered by excessive velocity' endif if (abs(ax)>accel_lim.or.abs(ay)>accel_lim) then dumpit=.true. write(stderr(),'("pe=",i3,x,a)') mpp_pe(),'Dump triggered by excessive acceleration' endif if (present(debug_flag)) then if (debug_flag) dumpit=.true. write(stderr(),'("pe=",i3,x,a)') mpp_pe(),'Debug dump flagged by arguments' endif if (dumpit) then 100 format('pe=',i3,a15,9(x,a8,es12.3)) 200 format('pe=',i3,a15,(x,a8,i12),9(x,a8,es12.3)) write(stderr(),200) mpp_pe(),'Starting pars:', & 'yr0=',berg%start_year, 'day0=',berg%start_day, & 'lon0=',berg%start_lon, 'lat0=',berg%start_lat, 'mass0=',berg%start_mass, & 'sclng=',berg%mass_scaling write(stderr(),100) mpp_pe(),'Geometry:', & 'M=',M, 'T=',T, 'D=',D, 'F=',F, 'W=',W, 'L=',L write(stderr(),100) mpp_pe(),'delta U:', & 'u(n)=',uvel0, 'u(*)=', uvel, 'u(n+1)=',uvel+dt*ax, 'del u=',dt*ax write(stderr(),100) mpp_pe(),'U terms', & 'f*v=',f_cori*vvel, & 'g*H_x=',-gravity*ssh_x, & 'wave*ua=',wave_rad*uwave, & 'd*(u-uo)=',-drag_ocn*(uvel-uo), & 'd*(u-ua)=',-drag_atm*(uvel-ua), & 'd*(u-ui)=',-drag_ice*(uvel-ui) write(stderr(),100) mpp_pe(),'U accel.', & 'axe=',axe, & 'ax=',ax, & 'ax(cori)=',detA*(A11*(f_cori*vvel)+A12*(-f_cori*uvel)), & 'ax(grav)=',detA*(A11*(-gravity*ssh_x)+A12*(-gravity*ssh_y)), & 'ax(wave)=',detA*(A11*(wave_rad*uwave)+A12*(wave_rad*vwave)), & 'ax(ocn)=',detA*(A11*(-drag_ocn*(uvel-uo))+A12*(-drag_ocn*(vvel-vo))), & 'ax(atm)=',detA*(A11*(-drag_atm*(uvel-ua))+A12*(-drag_atm*(vvel-va))), & 'ax(ice)=',detA*(A11*(-drag_ice*(uvel-ui))+A12*(-drag_ice*(vvel-vi))) write(stderr(),100) mpp_pe(),'delta V:', & 'v(n)=',vvel0, 'v(*)=', vvel, 'v(n+1)=',vvel+dt*ay, 'del v=',dt*ay write(stderr(),100) mpp_pe(),'V terms', & 'f*u=',-f_cori*uvel, & 'g*H_y=',-gravity*ssh_y, & 'wave*va=',wave_rad*vwave, & 'd*(v-vo)=',-drag_ocn*(vvel-vo), & 'd*(v-va)=',-drag_atm*(vvel-va), & 'd*(v-vi)=',-drag_ice*(vvel-vi) write(stderr(),100) mpp_pe(),'V accel. pe=', & 'aye=',aye, & 'ay=',ay, & 'ay(cori)=',detA*(-A12*(f_cori*vvel)+A11*(-f_cori*uvel)), & 'ay(grav)=',detA*(-A12*(-gravity*ssh_x)+A11*(-gravity*ssh_y)), & 'ay(wave)=',detA*(-A12*(wave_rad*uwave)+A11*(wave_rad*vwave)), & 'ay(ocn)=',detA*(-A12*(-drag_ocn*(uvel-uo))+A11*(-drag_ocn*(vvel-vo))), & 'ay(atm)=',detA*(-A12*(-drag_atm*(uvel-ua))+A11*(-drag_atm*(vvel-va))), & 'ay(ice)=',detA*(-A12*(-drag_ice*(uvel-ui))+A11*(-drag_ice*(vvel-vi))) write(stderr(),100) mpp_pe(),'Vel scales', & '|va-vo|=',sqrt((ua-uo)**2+(va-vo)**2), & '|vo-vb|=',sqrt((uvel-uo)**2+(vvel-vo)**2), & '|va-vb|=',sqrt((uvel-ua)**2+(vvel-va)**2), & '|vi-vb|=',sqrt((uvel-ui)**2+(vvel-vi)**2), & '|vb|=',sqrt((uvel)**2+(vvel)**2), & '|va|=',sqrt((ua)**2+(va)**2), & '|vo|=',sqrt((uo)**2+(vo)**2), & '|vi|=',sqrt((ui)**2+(vi)**2) write(stderr(),100) mpp_pe(),'Time scales', & 'f=',f_cori, 'wave_rad=',wave_rad, 'do=',drag_ocn, 'da=',drag_atm, 'di=',drag_ice write(stderr(),100) mpp_pe(),'u*', & 'd*=',lambda, & 'u*=',(drag_ocn*uo+drag_atm*ua+drag_ice*ui)/lambda, & 'uo*=',(drag_ocn*uo)/lambda, & 'ua*=',(drag_atm*ua)/lambda, & 'ui*=',(drag_ice*ui)/lambda write(stderr(),100) mpp_pe(),'v*', & 'd*=',lambda, & 'v*=',(drag_ocn*vo+drag_atm*va+drag_ice*vi)/lambda, & 'vo*=',(drag_ocn*vo)/lambda, & 'va*=',(drag_atm*va)/lambda, & 'vi*=',(drag_ice*vi)/lambda write(stderr(),100) mpp_pe(),'params', & 'a=',ampl, 'Lwl=',Lwavelength, 'Lcut=',Lcutoff, 'Ltop=',Ltop, 'hi=',hi, 'Cr=',Cr write(stderr(),100) mpp_pe(),'Position', & 'xi=',xi, 'yj=',yj, 'lat=',lat call dump_locfld(grd,i,j,grd%msk,'MSK') call dump_locfld(grd,i,j,grd%ssh,'SSH') call dump_locfld(grd,i,j,grd%sst,'SST') call dump_locvel(grd,i,j,grd%uo,'Uo') call dump_locvel(grd,i,j,grd%vo,'Vo') call dump_locvel(grd,i,j,grd%ua,'Ua') call dump_locvel(grd,i,j,grd%va,'Va') call dump_locvel(grd,i,j,grd%ui,'Ui') call dump_locvel(grd,i,j,grd%vi,'Vi') call dump_locfld(grd,i,j,grd%hi,'HI') call dump_locfld(grd,i,j,grd%cn,'CN') call dump_locvel(grd,i,j,grd%lon,'Lon') call dump_locvel(grd,i,j,grd%lat,'Lat') call print_berg(stderr(),berg,'diamonds, accel, large accel') endif contains subroutine dump_locfld(grd,i0,j0,A,lbl) ! Arguments type(icebergs_gridded), pointer :: grd integer, intent(in) :: i0, j0 real, dimension(grd%isd:grd%ied,grd%jsd:grd%jed), intent(in) :: A character(len=*) :: lbl ! Local variables integer :: i, j, ii, jj real :: B(-1:1,-1:1), fac do jj=-1,1 j=max(grd%jsd,min(grd%jed,jj+j0)) do ii=-1,1 i=max(grd%isd,min(grd%ied,ii+i0)) B(ii,jj)=A(i,j) if ((i.ne.ii+i0).or.(j.ne.jj+j0)) B(ii,jj)=-9.999999e-99 enddo enddo write(stderr(),'("pe=",i3,x,a8,3i12)') mpp_pe(),lbl,(i0+ii,ii=-1,1) do jj=1,-1,-1 write(stderr(),'("pe=",i3,x,i8,3es12.4)') mpp_pe(),j0+jj,(B(ii,jj),ii=-1,1) enddo end subroutine dump_locfld subroutine dump_locvel(grd,i0,j0,A,lbl) ! Arguments type(icebergs_gridded), pointer :: grd integer, intent(in) :: i0, j0 real, dimension(grd%isd:grd%ied,grd%jsd:grd%jed), intent(in) :: A character(len=*) :: lbl ! Local variables integer :: i, j, ii, jj real :: B(-1:0,-1:0), fac do jj=-1,0 j=max(grd%jsd,min(grd%jed,jj+j0)) do ii=-1,0 i=max(grd%isd,min(grd%ied,ii+i0)) B(ii,jj)=A(i,j) if ((i.ne.ii+i0).or.(j.ne.jj+j0)) B(ii,jj)=-9.999999e-99 enddo enddo write(stderr(),'("pe=",i3,x,a8,3i12)') mpp_pe(),lbl,(i0+ii,ii=-1,0) do jj=0,-1,-1 write(stderr(),'("pe=",i3,x,i8,3es12.4)') mpp_pe(),j0+jj,(B(ii,jj),ii=-1,0) enddo end subroutine dump_locvel end subroutine accel ! ############################################################################## subroutine thermodynamics(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables type(icebergs_gridded), pointer :: grd real :: M, T, W, L, SST, Vol, Ln, Wn, Tn, nVol, IC, Dn real :: Mv, Me, Mb, melt, dvo, dva, dM, Ss, dMe, dMb, dMv real :: Mnew, Mnew1, Mnew2 real :: Mbits, nMbits, dMbitsE, dMbitsM, Lbits, Abits, Mbb integer :: i,j type(iceberg), pointer :: this, next real, parameter :: perday=1./86400. ! For convenience grd=>bergs%grd this=>bergs%first do while(associated(this)) if (debug) call check_position(grd, this, 'thermodynamics (top)') call interp_flds(grd, this%ine, this%jne, this%xi, this%yj, this%uo, this%vo, & this%ui, this%vi, this%ua, this%va, this%ssh_x, this%ssh_y, this%sst, & this%cn, this%hi) SST=this%sst IC=min(1.,this%cn+bergs%sicn_shift) ! Shift sea-ice concentration M=this%mass T=this%thickness ! total thickness ! D=(bergs%rho_bergs/rho_seawater)*T ! draught (keel depth) ! F=T-D ! freeboard W=this%width L=this%length i=this%ine j=this%jne Vol=T*W*L ! Environment dvo=sqrt((this%uvel-this%uo)**2+(this%vvel-this%vo)**2) dva=sqrt((this%ua-this%uo)**2+(this%va-this%vo)**2) Ss=1.5*(dva**0.5)+0.1*dva ! Sea state ! Melt rates in m/s Mv=max( 7.62e-3*SST+1.29e-3*(SST**2), 0.) &! Buoyant convection at sides *perday ! convert to m/s Mb=max( 0.58*(dvo**0.8)*(SST+4.0)/(L**0.2), 0.) &! Basal turbulent melting *perday ! convert to m/s Me=max( 1./12.*(SST+2.)*Ss*(1+cos(pi*(IC**3))) ,0.) &! Wave erosion *perday ! convert to m/s if (bergs%use_operator_splitting) then ! Operator split update of volume/mass Tn=max(T-Mb*bergs%dt,0.) ! new total thickness (m) nVol=Tn*W*L ! new volume (m^3) Mnew1=(nVol/Vol)*M ! new mass (kg) dMb=M-Mnew1 ! mass lost to basal melting (>0) (kg) Ln=max(L-Mv*bergs%dt,0.) ! new length (m) Wn=max(W-Mv*bergs%dt,0.) ! new width (m) nVol=Tn*Wn*Ln ! new volume (m^3) Mnew2=(nVol/Vol)*M ! new mass (kg) dMv=Mnew1-Mnew2 ! mass lost to buoyant convection (>0) (kg) Ln=max(Ln-Me*bergs%dt,0.) ! new length (m) Wn=max(Wn-Me*bergs%dt,0.) ! new width (m) nVol=Tn*Wn*Ln ! new volume (m^3) Mnew=(nVol/Vol)*M ! new mass (kg) dMe=Mnew2-Mnew ! mass lost to erosion (>0) (kg) dM=M-Mnew ! mass lost to all erosion and melting (>0) (kg) else ! Update dimensions of berg Ln=max(L-(Mv+Me)*(bergs%dt),0.) ! (m) Wn=max(W-(Mv+Me)*(bergs%dt),0.) ! (m) Tn=max(T-Mb*(bergs%dt),0.) ! (m) ! Update volume and mass of berg nVol=Tn*Wn*Ln ! (m^3) Mnew=(nVol/Vol)*M ! (kg) dM=M-Mnew ! (kg) dMb=(M/Vol)*(W*L)*Mb*bergs%dt ! approx. mass loss to basal melting (kg) dMe=(M/Vol)*(T*(W+L))*Me*bergs%dt ! approx. mass lost to erosion (kg) dMv=(M/Vol)*(T*(W+L))*Mv*bergs%dt ! approx. mass loss to buoyant convection (kg) endif ! Bergy bits if (bergs%bergy_bit_erosion_fraction>0.) then Mbits=this%mass_of_bits ! mass of bergy bits (kg) dMbitsE=bergs%bergy_bit_erosion_fraction*dMe ! change in mass of bits (kg) nMbits=Mbits+dMbitsE ! add new bergy bits to mass (kg) Lbits=min(L,W,T,40.) ! assume bergy bits are smallest dimension or 40 meters Abits=(Mbits/bergs%rho_bergs)/Lbits ! Effective bottom area (assuming T=Lbits) Mbb=max( 0.58*(dvo**0.8)*(SST+2.0)/(Lbits**0.2), 0.) &! Basal turbulent melting (for bits) *perday ! convert to m/s Mbb=bergs%rho_bergs*Abits*Mbb ! in kg/s dMbitsM=min(Mbb*bergs%dt,nMbits) ! bergy bits mass lost to melting (kg) nMbits=nMbits-dMbitsM ! remove mass lost to bergy bits melt if (Mnew==0.) then ! if parent berg has completely melted then dMbitsM=dMbitsM+nMbits ! instantly melt all the bergy bits nMbits=0. endif else Abits=0. dMbitsE=0. dMbitsM=0. nMbits=this%mass_of_bits ! retain previous value incase non-zero endif ! Add melting to the grid and field diagnostics if (grd%area(i,j).ne.0.) then melt=(dM-(dMbitsE-dMbitsM))/bergs%dt ! kg/s grd%floating_melt(i,j)=grd%floating_melt(i,j)+melt/grd%area(i,j)*this%mass_scaling ! kg/m2/s melt=melt*this%heat_density ! kg/s x J/kg = J/s grd%calving_hflx(i,j)=grd%calving_hflx(i,j)+melt/grd%area(i,j)*this%mass_scaling ! W/m2 bergs%net_heat_to_ocean=bergs%net_heat_to_ocean+melt*this%mass_scaling*bergs%dt ! J melt=dM/bergs%dt ! kg/s grd%berg_melt(i,j)=grd%berg_melt(i,j)+melt/grd%area(i,j)*this%mass_scaling ! kg/m2/s melt=dMbitsE/bergs%dt ! mass flux into bergy bits in kg/s grd%bergy_src(i,j)=grd%bergy_src(i,j)+melt/grd%area(i,j)*this%mass_scaling ! kg/m2/s melt=dMbitsM/bergs%dt ! melt rate of bergy bits in kg/s grd%bergy_melt(i,j)=grd%bergy_melt(i,j)+melt/grd%area(i,j)*this%mass_scaling ! kg/m2/s if(grd%id_melt_buoy>0) then melt=dMb/bergs%dt ! melt rate due to buoyancy term in kg/s grd%melt_buoy(i,j)=grd%melt_buoy(i,j)+melt/grd%area(i,j)*this%mass_scaling ! kg/m2/s endif if(grd%id_melt_eros>0) then melt=dMe/bergs%dt ! erosion rate in kg/s grd%melt_eros(i,j)=grd%melt_eros(i,j)+melt/grd%area(i,j)*this%mass_scaling ! kg/m2/s endif if(grd%id_melt_conv>0) then melt=dMv/bergs%dt ! melt rate due to convection term in kg/s grd%melt_conv(i,j)=grd%melt_conv(i,j)+melt/grd%area(i,j)*this%mass_scaling ! kg/m2/s endif else write(stderr(),*) 'diamonds, thermodynamics: berg appears to have grounded!!!! PE=',mpp_pe(),i,j call print_berg(stderr(),this,'thermodynamics, grounded') if (associated(this%trajectory)) & write(stderr(),*) 'traj=',this%trajectory%lon,this%trajectory%lat write(stderr(),*) 'msk=',grd%msk(i,j),grd%area(i,j) call error_mesg('diamonds, thermodynamics', 'berg appears to have grounded!', FATAL) endif ! Rolling Dn=(bergs%rho_bergs/rho_seawater)*Tn ! draught (keel depth) if ( Dn>0. .and. max(Wn,Ln)W) this%mass=Mnew this%mass_of_bits=nMbits this%thickness=Tn this%width=min(Wn,Ln) this%length=max(Wn,Ln) next=>this%next ! Did berg completely melt? if (Mnew<=0.) then ! Delete the berg call move_trajectory(bergs, this) call delete_iceberg_from_list(bergs%first, this) bergs%nbergs_melted=bergs%nbergs_melted+1 else ! Diagnose mass distribution on grid if (grd%id_virtual_area>0) grd%virtual_area(i,j)=grd%virtual_area(i,j)+(Wn*Ln+Abits)*this%mass_scaling ! m^2 if (grd%id_mass>0 .or. bergs%add_weight_to_ocean) grd%mass(i,j)=grd%mass(i,j)+Mnew/grd%area(i,j)*this%mass_scaling ! kg/m2 if (grd%id_bergy_mass>0 .or. bergs%add_weight_to_ocean) grd%bergy_mass(i,j)=grd%bergy_mass(i,j)+nMbits/grd%area(i,j)*this%mass_scaling ! kg/m2 if (bergs%add_weight_to_ocean .and. .not. bergs%time_average_weight) & call spread_mass_across_ocean_cells(grd, i, j, this%xi, this%yj, Mnew, nMbits, this%mass_scaling) endif this=>next enddo end subroutine thermodynamics ! ############################################################################## subroutine spread_mass_across_ocean_cells(grd, i, j, x, y, Mberg, Mbits, scaling) ! Arguments type(icebergs_gridded), pointer :: grd integer, intent(in) :: i, j real, intent(in) :: x, y, Mberg, Mbits, scaling ! Local variables real :: xL, xC, xR, yD, yC, yU, Mass real :: yDxL, yDxC, yDxR, yCxL, yCxC, yCxR, yUxL, yUxC, yUxR Mass=(Mberg+Mbits)*scaling xL=min(0.5, max(0., 0.5-x)) xR=min(0.5, max(0., x-0.5)) xC=max(0., 1.-(xL+xR)) yD=min(0.5, max(0., 0.5-y)) yU=min(0.5, max(0., y-0.5)) yC=max(0., 1.-(yD+yU)) yDxL=yD*xL*grd%msk(i-1,j-1) yDxC=yD*xC*grd%msk(i ,j-1) yDxR=yD*xR*grd%msk(i+1,j-1) yCxL=yC*xL*grd%msk(i-1,j ) yCxR=yC*xR*grd%msk(i+1,j ) yUxL=yU*xL*grd%msk(i-1,j+1) yUxC=yU*xC*grd%msk(i ,j+1) yUxR=yU*xR*grd%msk(i+1,j+1) yCxC=1.-( ((yDxL+yUxR)+(yDxR+yUxL)) + ((yCxL+yCxR)+(yDxC+yUxC)) ) grd%mass_on_ocean(i,j,1)=grd%mass_on_ocean(i,j,1)+yDxL*Mass grd%mass_on_ocean(i,j,2)=grd%mass_on_ocean(i,j,2)+yDxC*Mass grd%mass_on_ocean(i,j,3)=grd%mass_on_ocean(i,j,3)+yDxR*Mass grd%mass_on_ocean(i,j,4)=grd%mass_on_ocean(i,j,4)+yCxL*Mass grd%mass_on_ocean(i,j,5)=grd%mass_on_ocean(i,j,5)+yCxC*Mass grd%mass_on_ocean(i,j,6)=grd%mass_on_ocean(i,j,6)+yCxR*Mass grd%mass_on_ocean(i,j,7)=grd%mass_on_ocean(i,j,7)+yUxL*Mass grd%mass_on_ocean(i,j,8)=grd%mass_on_ocean(i,j,8)+yUxC*Mass grd%mass_on_ocean(i,j,9)=grd%mass_on_ocean(i,j,9)+yUxR*Mass end subroutine spread_mass_across_ocean_cells ! ############################################################################## subroutine interp_flds(grd, i, j, xi, yj, uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi) ! Arguments type(icebergs_gridded), pointer :: grd integer, intent(in) :: i, j real, intent(in) :: xi, yj real, intent(out) :: uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi ! Local variables real :: cos_rot, sin_rot #ifdef USE_OLD_SSH_GRADIENT real :: dxm, dx0, dxp #endif real :: hxm, hxp real, parameter :: ssh_coast=0.00 cos_rot=bilin(grd, grd%cos, i, j, xi, yj) sin_rot=bilin(grd, grd%sin, i, j, xi, yj) uo=bilin(grd, grd%uo, i, j, xi, yj) vo=bilin(grd, grd%vo, i, j, xi, yj) ui=bilin(grd, grd%ui, i, j, xi, yj) vi=bilin(grd, grd%vi, i, j, xi, yj) ua=bilin(grd, grd%ua, i, j, xi, yj) va=bilin(grd, grd%va, i, j, xi, yj) ! These fields are cell centered (A-grid) and would ! best be interpolated using PLM. For now we use PCM! sst=grd%sst(i,j) ! A-grid cn=grd%cn(i,j) ! A-grid hi=grd%hi(i,j) ! A-grid ! Estimate SSH gradient in X direction #ifdef USE_OLD_SSH_GRADIENT dxp=0.5*(grd%dx(i+1,j)+grd%dx(i+1,j-1)) dx0=0.5*(grd%dx(i,j)+grd%dx(i,j-1)) dxm=0.5*(grd%dx(i-1,j)+grd%dx(i-1,j-1)) hxm=2.*(grd%ssh(i,j)-grd%ssh(i-1,j))/(dx0+dxm)*grd%msk(i-1,j) & +(-ssh_coast)/(dx0+dxm)*(1.-grd%msk(i-1,j)) ! force to drive bergs away from coasts hxp=2.*(grd%ssh(i+1,j)-grd%ssh(i,j))/(dx0+dxp)*grd%msk(i+1,j) & +(+ssh_coast)/(dx0+dxp)*(1.-grd%msk(i+1,j)) ! force to drive bergs away from coasts #else if (yj>=0.5) then hxp=(yj-0.5)*ddx_ssh(grd,i ,j+1)+(1.5-yj)*ddx_ssh(grd,i ,j ) hxm=(yj-0.5)*ddx_ssh(grd,i-1,j+1)+(1.5-yj)*ddx_ssh(grd,i-1,j ) else hxp=(yj+0.5)*ddx_ssh(grd,i ,j )+(0.5-yj)*ddx_ssh(grd,i ,j-1) hxm=(yj+0.5)*ddx_ssh(grd,i-1,j )+(0.5-yj)*ddx_ssh(grd,i-1,j-1) endif #endif ! ssh_x is at the u-point on a C-grid ssh_x=xi*hxp+(1.-xi)*hxm ! Estimate SSH gradient in Y direction #ifdef USE_OLD_SSH_GRADIENT dxp=0.5*(grd%dy(i,j+1)+grd%dy(i-1,j+1)) dx0=0.5*(grd%dy(i,j)+grd%dy(i-1,j)) dxm=0.5*(grd%dy(i,j-1)+grd%dy(i-1,j-1)) hxm=2.*(grd%ssh(i,j)-grd%ssh(i,j-1))/(dx0+dxm)*grd%msk(i,j-1) & +(-ssh_coast)/(dx0+dxm)*(1.-grd%msk(i,j-1)) ! force to drive bergs away from coasts hxp=2.*(grd%ssh(i,j+1)-grd%ssh(i,j))/(dx0+dxp)*grd%msk(i,j+1) & +(+ssh_coast)/(dx0+dxp)*(1.-grd%msk(i,j+1)) ! force to drive bergs away from coasts #else if (xi>=0.5) then hxp=(xi-0.5)*ddy_ssh(grd,i+1,j )+(1.5-xi)*ddy_ssh(grd,i ,j ) hxm=(xi-0.5)*ddy_ssh(grd,i+1,j-1)+(1.5-xi)*ddy_ssh(grd,i ,j-1) else hxp=(xi+0.5)*ddy_ssh(grd,i ,j )+(0.5-xi)*ddy_ssh(grd,i-1,j ) hxm=(xi+0.5)*ddy_ssh(grd,i ,j-1)+(0.5-xi)*ddy_ssh(grd,i-1,j-1) endif #endif ! ssh_y is at the v-point on a C-grid ssh_y=yj*hxp+(1.-yj)*hxm ! Rotate vectors from local grid to lat/lon coordinates call rotate(uo, vo, cos_rot, sin_rot) call rotate(ui, vi, cos_rot, sin_rot) call rotate(ua, va, cos_rot, sin_rot) call rotate(ssh_x, ssh_y, cos_rot, sin_rot) contains real function ddx_ssh(grd,i,j) ! Arguments type(icebergs_gridded), pointer :: grd integer, intent(in) :: i, j ! Local variables real :: dxp,dx0 dxp=0.5*(grd%dx(i+1,j)+grd%dx(i+1,j-1)) dx0=0.5*(grd%dx(i,j)+grd%dx(i,j-1)) ddx_ssh=2.*(grd%ssh(i+1,j)-grd%ssh(i,j))/(dx0+dxp)*grd%msk(i+1,j)*grd%msk(i,j) end function ddx_ssh real function ddy_ssh(grd,i,j) ! Arguments type(icebergs_gridded), pointer :: grd integer, intent(in) :: i, j ! Local variables real :: dyp,dy0 dyp=0.5*(grd%dy(i,j+1)+grd%dy(i-1,j+1)) dy0=0.5*(grd%dy(i,j)+grd%dy(i-1,j)) ddy_ssh=2.*(grd%ssh(i,j+1)-grd%ssh(i,j))/(dy0+dyp)*grd%msk(i,j+1)*grd%msk(i,j) end function ddy_ssh subroutine rotate(u, v, cos_rot, sin_rot) ! Arguments real, intent(inout) :: u, v real, intent(in) :: cos_rot, sin_rot ! Local variables real :: u_old, v_old u_old=u v_old=v u=cos_rot*u_old+sin_rot*v_old v=cos_rot*v_old-sin_rot*u_old end subroutine rotate end subroutine interp_flds ! ############################################################################## real function bilin(grd, fld, i, j, xi, yj) ! Arguments type(icebergs_gridded), pointer :: grd real, intent(in) :: fld(grd%isd:grd%ied,grd%jsd:grd%jed), xi, yj integer, intent(in) :: i, j ! Local variables bilin=(fld(i,j )*(1.-xi)+fld(i-1,j )*xi)*(1.-yj) & +(fld(i,j-1)*(1.-xi)+fld(i-1,j-1)*xi)*yj end function bilin ! ############################################################################## subroutine icebergs_run(bergs, time, calving, uo, vo, ui, vi, tauxa, tauya, ssh, sst, calving_hflx, cn, hi) ! Arguments type(icebergs), pointer :: bergs type(time_type), intent(in) :: time real, dimension(:,:), intent(inout) :: calving, calving_hflx real, dimension(:,:), intent(in) :: uo, vo, ui, vi, tauxa, tauya, ssh, sst, cn, hi ! Local variables integer :: iyr, imon, iday, ihr, imin, isec, k type(icebergs_gridded), pointer :: grd logical :: lerr, sample_traj, lbudget, lverbose real :: unused_calving, tmpsum, grdd_berg_mass, grdd_bergy_mass call mpp_clock_begin(bergs%clock) call mpp_clock_begin(bergs%clock_int) ! For convenience grd=>bergs%grd grd%floating_melt(:,:)=0. grd%berg_melt(:,:)=0. grd%melt_buoy(:,:)=0. grd%melt_eros(:,:)=0. grd%melt_conv(:,:)=0. grd%bergy_src(:,:)=0. grd%bergy_melt(:,:)=0. grd%bergy_mass(:,:)=0. grd%mass(:,:)=0. if (bergs%add_weight_to_ocean) grd%mass_on_ocean(:,:,:)=0. grd%virtual_area(:,:)=0. ! Manage time call get_date(time, iyr, imon, iday, ihr, imin, isec) bergs%current_year=iyr bergs%current_yearday=float(imon-1)*31.+float(iday-1)+(float(ihr)+(float(imin)+float(isec)/60.)/60.)/24. ! Turn on sampling of trajectories, verbosity, budgets sample_traj=.false. if (bergs%traj_sample_hrs>0 .and. mod(24*iday+ihr,bergs%traj_sample_hrs).eq.0) sample_traj=.true. lverbose=.false. if (bergs%verbose_hrs>0 .and. mod(24*iday+ihr,bergs%verbose_hrs).eq.0) lverbose=verbose lbudget=.false. if (bergs%verbose_hrs>0 .and. mod(24*iday+ihr,bergs%verbose_hrs).eq.0) lbudget=budget if (mpp_pe()==mpp_root_pe().and.lverbose) write(*,'(a,3i5,a,3i5,a,i5,f8.3)') & 'diamonds: y,m,d=',iyr, imon, iday,' h,m,s=', ihr, imin, isec, & ' yr,yrdy=', bergs%current_year, bergs%current_yearday ! Adapt calving flux from coupler for use here !call sanitize_field(grd%calving,1.e20) tmpsum=sum( calving(:,:)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%net_calving_received=bergs%net_calving_received+tmpsum*bergs%dt grd%calving(grd%isc:grd%iec,grd%jsc:grd%jec)=calving(:,:) ! Units of kg/m2/s grd%calving(:,:)=grd%calving(:,:)*grd%msk(:,:)*grd%area(:,:) ! Convert to kg/s tmpsum=sum( grd%calving(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%net_incoming_calving=bergs%net_incoming_calving+tmpsum*bergs%dt if (grd%id_calving>0) & lerr=send_data(grd%id_calving, grd%calving(grd%isc:grd%iec,grd%jsc:grd%jec), Time) ! Adapt calving heat flux from coupler grd%calving_hflx(grd%isc:grd%iec,grd%jsc:grd%jec)=calving_hflx(:,:) ! Units of W/m2 grd%calving_hflx(:,:)=grd%calving_hflx(:,:)*grd%msk(:,:) ! Mask (just in case) if (grd%id_calving_hflx_in>0) & lerr=send_data(grd%id_calving_hflx_in, grd%calving_hflx(grd%isc:grd%iec,grd%jsc:grd%jec), Time) tmpsum=sum( grd%calving_hflx(grd%isc:grd%iec,grd%jsc:grd%jec)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%net_incoming_calving_heat=bergs%net_incoming_calving_heat+tmpsum*bergs%dt ! Units of J ! Copy ocean flow (resides on B grid) grd%uo(grd%isc-1:grd%iec+1,grd%jsc-1:grd%jec+1)=uo(:,:) grd%vo(grd%isc-1:grd%iec+1,grd%jsc-1:grd%jec+1)=vo(:,:) call mpp_update_domains(grd%uo, grd%vo, grd%domain, gridtype=BGRID_NE) ! Copy ice flow (resides on B grid) grd%ui(grd%isc-1:grd%iec+1,grd%jsc-1:grd%jec+1)=ui(:,:) grd%vi(grd%isc-1:grd%iec+1,grd%jsc-1:grd%jec+1)=vi(:,:) call mpp_update_domains(grd%ui, grd%vi, grd%domain, gridtype=BGRID_NE) ! Copy atmospheric stress (resides on A grid) grd%ua(grd%isc:grd%iec,grd%jsc:grd%jec)=tauxa(:,:) ! Note rough conversion from stress to speed grd%va(grd%isc:grd%iec,grd%jsc:grd%jec)=tauya(:,:) ! Note rough conversion from stress to speed call invert_tau_for_du(grd%ua, grd%va) ! Note rough conversion from stress to speed !grd%ua(grd%isc:grd%iec,grd%jsc:grd%jec)=sign(sqrt(abs(tauxa(:,:))/0.01),tauxa(:,:)) ! Note rough conversion from stress to speed !grd%va(grd%isc:grd%iec,grd%jsc:grd%jec)=sign(sqrt(abs(tauya(:,:))/0.01),tauya(:,:)) ! Note rough conversion from stress to speed call mpp_update_domains(grd%ua, grd%va, grd%domain, gridtype=BGRID_NE) ! Copy sea surface height and temperature(resides on A grid) grd%ssh(grd%isc-1:grd%iec+1,grd%jsc-1:grd%jec+1)=ssh(:,:) call mpp_update_domains(grd%ssh, grd%domain) grd%sst(grd%isc:grd%iec,grd%jsc:grd%jec)=sst(:,:)-273.15 ! Note convert from Kelvin to Celsius call mpp_update_domains(grd%sst, grd%domain) ! Copy sea-ice concentration and thickness (resides on A grid) grd%cn(grd%isc-1:grd%iec+1,grd%jsc-1:grd%jec+1)=cn(:,:) call mpp_update_domains(grd%cn, grd%domain) grd%hi(grd%isc-1:grd%iec+1,grd%jsc-1:grd%jec+1)=hi(:,:) call mpp_update_domains(grd%hi, grd%domain) if (debug) call bergs_chksum(bergs, 'run bergs (top)') if (debug) call checksum_gridded(bergs%grd, 'top of s/r run') ! Accumulate ice from calving call accumulate_calving(bergs) if (grd%id_accum>0) then grd%tmp(grd%isc:grd%iec,grd%jsc:grd%jec)=calving(:,:) grd%tmp(:,:)=grd%tmp(:,:)*grd%msk(:,:)*grd%area(:,:)-grd%calving(:,:) lerr=send_data(grd%id_accum, grd%tmp(grd%isc:grd%iec,grd%jsc:grd%jec), Time) endif if (grd%id_unused>0) & lerr=send_data(grd%id_unused, grd%calving(grd%isc:grd%iec,grd%jsc:grd%jec), Time) unused_calving=sum( grd%calving(grd%isc:grd%iec,grd%jsc:grd%jec) ) call mpp_clock_end(bergs%clock_int) call mpp_clock_begin(bergs%clock_cal) ! Calve excess stored ice into icebergs call calve_icebergs(bergs) if (debug) call bergs_chksum(bergs, 'run bergs (calved)') if (debug) call checksum_gridded(bergs%grd, 's/r run after calving') call mpp_clock_end(bergs%clock_cal) ! For each berg, evolve call mpp_clock_begin(bergs%clock_mom) if (associated(bergs%first)) call evolve_icebergs(bergs) if (debug) call bergs_chksum(bergs, 'run bergs (evolved)',ignore_halo_violation=.true.) if (debug) call checksum_gridded(bergs%grd, 's/r run after evolve') call mpp_clock_end(bergs%clock_mom) ! Send bergs to other PEs call mpp_clock_begin(bergs%clock_com) call send_bergs_to_other_pes(bergs) if (debug) call bergs_chksum(bergs, 'run bergs (exchanged)') if (debug) call checksum_gridded(bergs%grd, 's/r run after exchange') call mpp_clock_end(bergs%clock_com) ! Ice berg thermodynamics (melting) + rolling call mpp_clock_begin(bergs%clock_the) if (associated(bergs%first)) call thermodynamics(bergs) if (debug) call bergs_chksum(bergs, 'run bergs (thermo)') if (debug) call checksum_gridded(bergs%grd, 's/r run after thermodynamics') call mpp_clock_end(bergs%clock_the) ! For each berg, record call mpp_clock_begin(bergs%clock_dia) if (sample_traj.and.associated(bergs%first)) call record_posn(bergs) ! Gridded diagnostics if (grd%id_uo>0) & lerr=send_data(grd%id_uo, grd%uo(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_vo>0) & lerr=send_data(grd%id_vo, grd%vo(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_ui>0) & lerr=send_data(grd%id_ui, grd%ui(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_vi>0) & lerr=send_data(grd%id_vi, grd%vi(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_ua>0) & lerr=send_data(grd%id_ua, grd%ua(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_va>0) & lerr=send_data(grd%id_va, grd%va(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_sst>0) & lerr=send_data(grd%id_sst, grd%sst(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_cn>0) & lerr=send_data(grd%id_cn, grd%cn(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_hi>0) & lerr=send_data(grd%id_hi, grd%hi(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_floating_melt>0) & lerr=send_data(grd%id_floating_melt, grd%floating_melt(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_berg_melt>0) & lerr=send_data(grd%id_berg_melt, grd%berg_melt(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_melt_buoy>0) & lerr=send_data(grd%id_melt_buoy, grd%melt_buoy(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_melt_eros>0) & lerr=send_data(grd%id_melt_eros, grd%melt_eros(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_melt_conv>0) & lerr=send_data(grd%id_melt_conv, grd%melt_conv(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_virtual_area>0) & lerr=send_data(grd%id_virtual_area, grd%virtual_area(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_bergy_src>0) & lerr=send_data(grd%id_bergy_src, grd%bergy_src(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_bergy_melt>0) & lerr=send_data(grd%id_bergy_melt, grd%bergy_melt(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_bergy_mass>0) & lerr=send_data(grd%id_bergy_mass, grd%bergy_mass(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_mass>0) & lerr=send_data(grd%id_mass, grd%mass(grd%isc:grd%iec,grd%jsc:grd%jec), Time) if (grd%id_stored_ice>0) & lerr=send_data(grd%id_stored_ice, grd%stored_ice(grd%isc:grd%iec,grd%jsc:grd%jec,:), Time) if (grd%id_real_calving>0) & lerr=send_data(grd%id_real_calving, grd%real_calving(grd%isc:grd%iec,grd%jsc:grd%jec,:), Time) ! Dump icebergs to screen if (really_debug) call print_bergs(stderr(),bergs,'icebergs_run, status') call mpp_clock_end(bergs%clock_dia) ! Return what ever calving we did not use and additional icebergs melt call mpp_clock_begin(bergs%clock_int) if (.not. bergs%passive_mode) then where (grd%area(grd%isc:grd%iec,grd%jsc:grd%jec)>0.) calving(:,:)=grd%calving(grd%isc:grd%iec,grd%jsc:grd%jec)/grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) & +grd%floating_melt(grd%isc:grd%iec,grd%jsc:grd%jec) elsewhere calving(:,:)=0. end where calving_hflx(:,:)=grd%calving_hflx(grd%isc:grd%iec,grd%jsc:grd%jec) endif call mpp_clock_end(bergs%clock_int) ! Diagnose budgets call mpp_clock_begin(bergs%clock_dia) tmpsum=sum( grd%floating_melt(grd%isc:grd%iec,grd%jsc:grd%jec)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%net_melt=bergs%net_melt+tmpsum*bergs%dt bergs%net_outgoing_calving=bergs%net_outgoing_calving+(unused_calving+tmpsum)*bergs%dt tmpsum=sum( grd%berg_melt(grd%isc:grd%iec,grd%jsc:grd%jec)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%berg_melt=bergs%berg_melt+tmpsum*bergs%dt tmpsum=sum( grd%bergy_src(grd%isc:grd%iec,grd%jsc:grd%jec)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%bergy_src=bergs%bergy_src+tmpsum*bergs%dt tmpsum=sum( grd%bergy_melt(grd%isc:grd%iec,grd%jsc:grd%jec)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%bergy_melt=bergs%bergy_melt+tmpsum*bergs%dt tmpsum=sum( calving(:,:)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%net_calving_returned=bergs%net_calving_returned+tmpsum*bergs%dt tmpsum=sum( grd%calving_hflx(grd%isc:grd%iec,grd%jsc:grd%jec)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%net_outgoing_calving_heat=bergs%net_outgoing_calving_heat+tmpsum*bergs%dt ! Units of J if (lbudget) then bergs%stored_end=sum( grd%stored_ice(grd%isc:grd%iec,grd%jsc:grd%jec,:) ) bergs%stored_heat_end=sum( grd%stored_heat(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%floating_mass_end=sum_mass(bergs%first) bergs%icebergs_mass_end=sum_mass(bergs%first,justbergs=.true.) bergs%bergy_mass_end=sum_mass(bergs%first,justbits=.true.) bergs%floating_heat_end=sum_heat(bergs%first) grd%tmpc(:,:)=0.; call mpp_clock_end(bergs%clock); call mpp_clock_end(bergs%clock_dia) ! To enable calling of public s/r call icebergs_incr_mass(bergs, grd%tmpc) call mpp_clock_begin(bergs%clock_dia); call mpp_clock_begin(bergs%clock) ! To enable calling of public s/r bergs%returned_mass_on_ocean=sum( grd%tmpc(grd%isc:grd%iec,grd%jsc:grd%jec)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ) bergs%nbergs_end=count_bergs(bergs) call mpp_sum(bergs%stored_end) call mpp_sum(bergs%stored_heat_end) call mpp_sum(bergs%floating_mass_end) call mpp_sum(bergs%icebergs_mass_end) call mpp_sum(bergs%bergy_mass_end) call mpp_sum(bergs%floating_heat_end) call mpp_sum(bergs%returned_mass_on_ocean) call mpp_sum(bergs%nbergs_end) call mpp_sum(bergs%nbergs_calved) do k=1,nclasses; call mpp_sum(bergs%nbergs_calved_by_class(k)); enddo call mpp_sum(bergs%nbergs_melted) call mpp_sum(bergs%nspeeding_tickets) call mpp_sum(bergs%net_calving_returned) call mpp_sum(bergs%net_outgoing_calving) call mpp_sum(bergs%net_calving_received) call mpp_sum(bergs%net_incoming_calving) call mpp_sum(bergs%net_incoming_calving_heat) call mpp_sum(bergs%net_incoming_calving_heat_used) call mpp_sum(bergs%net_outgoing_calving_heat) call mpp_sum(bergs%net_calving_used) call mpp_sum(bergs%net_calving_to_bergs) call mpp_sum(bergs%net_heat_to_bergs) call mpp_sum(bergs%net_heat_to_ocean) call mpp_sum(bergs%net_melt) call mpp_sum(bergs%berg_melt) call mpp_sum(bergs%bergy_src) call mpp_sum(bergs%bergy_melt) grdd_berg_mass=sum( grd%mass(grd%isc:grd%iec,grd%jsc:grd%jec)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ); call mpp_sum(grdd_berg_mass) grdd_bergy_mass=sum( grd%bergy_mass(grd%isc:grd%iec,grd%jsc:grd%jec)*grd%area(grd%isc:grd%iec,grd%jsc:grd%jec) ); call mpp_sum(grdd_bergy_mass) if (mpp_pe().eq.mpp_root_pe()) then 100 format("diamonds: ",a19,3(a18,"=",es14.7,x,a2,:,","),a12,i8) 200 format("diamonds: ",a19,10(a18,"=",es14.7,x,a2,:,",")) call report_state('stored ice','kg','',bergs%stored_start,'',bergs%stored_end,'') call report_state('floating','kg','',bergs%floating_mass_start,'',bergs%floating_mass_end,'',bergs%nbergs_end) call report_state('icebergs','kg','',bergs%icebergs_mass_start,'',bergs%icebergs_mass_end,'') call report_state('bits','kg','',bergs%bergy_mass_start,'',bergs%bergy_mass_end,'') call report_istate('berg #','',bergs%nbergs_start,'',bergs%nbergs_end,'') call report_ibudget('berg #','calved',bergs%nbergs_calved, & 'melted',bergs%nbergs_melted, & '#',bergs%nbergs_start,bergs%nbergs_end) call report_budget('stored mass','kg','calving used',bergs%net_calving_used, & 'bergs',bergs%net_calving_to_bergs, & 'stored mass',bergs%stored_start,bergs%stored_end) call report_budget('floating mass','kg','calving used',bergs%net_calving_to_bergs, & 'bergs',bergs%net_melt, & 'stored mass',bergs%floating_mass_start,bergs%floating_mass_end) call report_budget('berg mass','kg','calving',bergs%net_calving_to_bergs, & 'melt+eros',bergs%berg_melt, & 'berg mass',bergs%icebergs_mass_start,bergs%icebergs_mass_end) call report_budget('bits mass','kg','eros used',bergs%bergy_src, & 'bergs',bergs%bergy_melt, & 'stored mass',bergs%bergy_mass_start,bergs%bergy_mass_end) call report_budget('net mass','kg','recvd',bergs%net_calving_received, & 'rtrnd',bergs%net_calving_returned, & 'net mass',bergs%stored_start+bergs%floating_mass_start, & bergs%stored_end+bergs%floating_mass_end) call report_consistant('iceberg mass','kg','gridded',grdd_berg_mass,'bergs',bergs%icebergs_mass_end) call report_consistant('bits mass','kg','gridded',grdd_bergy_mass,'bits',bergs%bergy_mass_end) call report_consistant('wieght','kg','returned',bergs%returned_mass_on_ocean,'floating',bergs%floating_mass_end) call report_state('net heat','J','',bergs%stored_heat_start+bergs%floating_heat_start,'',bergs%stored_heat_end+bergs%floating_heat_end,'') call report_state('stored heat','J','',bergs%stored_heat_start,'',bergs%stored_heat_end,'') call report_state('floating heat','J','',bergs%floating_heat_start,'',bergs%floating_heat_end,'') call report_budget('net heat','J','net heat',bergs%net_incoming_calving_heat, & 'net heat',bergs%net_outgoing_calving_heat, & 'net heat',bergs%stored_heat_start+bergs%floating_heat_start, & bergs%stored_heat_end+bergs%floating_heat_end) call report_budget('stored heat','J','calving used',bergs%net_incoming_calving_heat_used, & 'bergs',bergs%net_heat_to_bergs, & 'net heat',bergs%stored_heat_start,bergs%stored_heat_end) call report_budget('flting heat','J','calved',bergs%net_heat_to_bergs, & 'melt',bergs%net_heat_to_ocean, & 'net heat',bergs%floating_heat_start,bergs%floating_heat_end) if (debug) then call report_consistant('top interface','kg','from SIS',bergs%net_incoming_calving,'seen by diamonds',bergs%net_calving_received) call report_consistant('bot interface','kg','sent',bergs%net_outgoing_calving,'seen by SIS',bergs%net_calving_returned) endif write(*,'("diamonds: calved by class = ",i,20(",",i))') (bergs%nbergs_calved_by_class(k),k=1,nclasses) if (bergs%nspeeding_tickets>0) write(*,'("diamonds: speeding tickets issued = ",i)') bergs%nspeeding_tickets endif bergs%nbergs_start=bergs%nbergs_end bergs%stored_start=bergs%stored_end bergs%nbergs_melted=0 bergs%nbergs_calved=0 bergs%nbergs_calved_by_class(:)=0 bergs%nspeeding_tickets=0 bergs%stored_heat_start=bergs%stored_heat_end bergs%floating_heat_start=bergs%floating_heat_end bergs%floating_mass_start=bergs%floating_mass_end bergs%icebergs_mass_start=bergs%icebergs_mass_end bergs%bergy_mass_start=bergs%bergy_mass_end bergs%net_calving_used=0. bergs%net_calving_to_bergs=0. bergs%net_heat_to_bergs=0. bergs%net_heat_to_ocean=0. bergs%net_calving_received=0. bergs%net_calving_returned=0. bergs%net_incoming_calving=0. bergs%net_outgoing_calving=0. bergs%net_incoming_calving_heat=0. bergs%net_incoming_calving_heat_used=0. bergs%net_outgoing_calving_heat=0. bergs%net_melt=0. bergs%berg_melt=0. bergs%bergy_melt=0. bergs%bergy_src=0. endif if (debug) call bergs_chksum(bergs, 'run bergs (bot)') if (debug) call checksum_gridded(bergs%grd, 'end of s/r run') call mpp_clock_end(bergs%clock_dia) call mpp_clock_end(bergs%clock) contains subroutine report_state(budgetstr,budgetunits,startstr,startval,endstr,endval,delstr,nbergs) ! Arguments character*(*), intent(in) :: budgetstr, budgetunits, startstr, endstr, delstr real, intent(in) :: startval, endval integer, intent(in), optional :: nbergs ! Local variables if (present(nbergs)) then write(*,100) budgetstr//' state:', & startstr//' start',startval,budgetunits, & endstr//' end',endval,budgetunits, & 'Delta '//delstr,endval-startval,budgetunits, & '# of bergs',nbergs else write(*,100) budgetstr//' state:', & startstr//' start',startval,budgetunits, & endstr//' end',endval,budgetunits, & delstr//'Delta',endval-startval,budgetunits endif 100 format("diamonds: ",a19,3(a18,"=",es14.7,x,a2,:,","),a12,i8) end subroutine report_state subroutine report_consistant(budgetstr,budgetunits,startstr,startval,endstr,endval) ! Arguments character*(*), intent(in) :: budgetstr, budgetunits, startstr, endstr real, intent(in) :: startval, endval ! Local variables write(*,200) budgetstr//' check:', & startstr,startval,budgetunits, & endstr,endval,budgetunits, & 'error',(endval-startval)/((endval+startval)+1e-30),'nd' 200 format("diamonds: ",a19,10(a18,"=",es14.7,x,a2,:,",")) end subroutine report_consistant subroutine report_budget(budgetstr,budgetunits,instr,inval,outstr,outval,delstr,startval,endval) ! Arguments character*(*), intent(in) :: budgetstr, budgetunits, instr, outstr, delstr real, intent(in) :: inval, outval, startval, endval ! Local variables write(*,200) budgetstr//' budget:', & instr//' in',inval,budgetunits, & outstr//' out',outval,budgetunits, & 'Delta '//delstr,inval-outval,budgetunits, & 'error',((endval-startval)-(inval-outval))/max(1.e-30,max(abs(endval-startval),abs(inval-outval))),'nd' 200 format("diamonds: ",a19,10(a18,"=",es14.7,x,a2,:,",")) end subroutine report_budget subroutine report_istate(budgetstr,startstr,startval,endstr,endval,delstr) ! Arguments character*(*), intent(in) :: budgetstr, startstr, endstr, delstr integer, intent(in) :: startval, endval ! Local variables write(*,100) budgetstr//' state:', & startstr//' start',startval, & endstr//' end',endval, & delstr//'Delta',endval-startval 100 format("diamonds: ",a19,3(a18,"=",i14,x,:,",")) end subroutine report_istate subroutine report_ibudget(budgetstr,instr,inval,outstr,outval,delstr,startval,endval) ! Arguments character*(*), intent(in) :: budgetstr, instr, outstr, delstr integer, intent(in) :: inval, outval, startval, endval ! Local variables write(*,200) budgetstr//' budget:', & instr//' in',inval, & outstr//' out',outval, & 'Delta '//delstr,inval-outval, & 'error',((endval-startval)-(inval-outval)) 200 format("diamonds: ",a19,10(a18,"=",i14,x,:,",")) end subroutine report_ibudget end subroutine icebergs_run ! ############################################################################## subroutine icebergs_incr_mass(bergs, mass) ! Arguments type(icebergs), pointer :: bergs real, dimension(bergs%grd%isc:bergs%grd%iec,bergs%grd%jsc:bergs%grd%jec), intent(inout) :: mass ! Local variables integer :: i, j type(icebergs_gridded), pointer :: grd real :: dmda if (.not. associated(bergs)) return if (.not. bergs%add_weight_to_ocean) return call mpp_clock_begin(bergs%clock) call mpp_clock_begin(bergs%clock_int) ! For convenience grd=>bergs%grd ! Add iceberg+bits mass field to non-haloed SIS field (kg/m^2) !mass(:,:)=mass(:,:)+( grd%mass(grd%isc:grd%iec,grd%jsc:grd%jec) & ! + grd%bergy_mass(grd%isc:grd%iec,grd%jsc:grd%jec) ) if (debug) then bergs%grd%tmp(:,:)=0.; bergs%grd%tmp(grd%isc:grd%iec,grd%jsc:grd%jec)=mass call grd_chksum2(bergs%grd, bergs%grd%tmp, 'mass in (incr)') endif call mpp_update_domains(grd%mass_on_ocean, grd%domain) do j=grd%jsc, grd%jec; do i=grd%isc, grd%iec dmda=grd%mass_on_ocean(i,j,5) & + ( ( (grd%mass_on_ocean(i-1,j-1,9)+grd%mass_on_ocean(i+1,j+1,1)) & + (grd%mass_on_ocean(i+1,j-1,7)+grd%mass_on_ocean(i-1,j+1,3)) ) & + ( (grd%mass_on_ocean(i-1,j ,6)+grd%mass_on_ocean(i+1,j ,4)) & + (grd%mass_on_ocean(i ,j-1,8)+grd%mass_on_ocean(i ,j+1,2)) ) ) if (grd%area(i,j)>0) dmda=dmda/grd%area(i,j)*grd%msk(i,j) if (.not. bergs%passive_mode) mass(i,j)=mass(i,j)+dmda enddo; enddo if (debug) then call grd_chksum3(bergs%grd, bergs%grd%mass_on_ocean, 'mass bergs (incr)') bergs%grd%tmp(:,:)=0.; bergs%grd%tmp(grd%isc:grd%iec,grd%jsc:grd%jec)=mass call grd_chksum2(bergs%grd, bergs%grd%tmp, 'mass out (incr)') endif call mpp_clock_end(bergs%clock_int) call mpp_clock_end(bergs%clock) end subroutine icebergs_incr_mass ! ############################################################################## subroutine accumulate_calving(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables type(icebergs_gridded), pointer :: grd real :: remaining_dist, net_calving_used integer :: k, i, j logical, save :: first_call=.true. ! For convenience grd=>bergs%grd ! This is a hack to simplify initialization if (first_call.and..not.bergs%restarted) then first_call=.false. !do k=1, nclasses ! where (grd%calving==0.) grd%stored_ice(:,:,k)=0. !enddo bergs%stored_start=sum( grd%stored_ice(grd%isc:grd%iec,grd%jsc:grd%jec,:) ) call mpp_sum( bergs%stored_start ) if (mpp_pe().eq.mpp_root_pe()) write(*,'(a,es13.6,a)') & 'diamonds, accumulate_calving: initial stored mass=',bergs%stored_start,' kg' do j=grd%jsc,grd%jec; do i=grd%isc,grd%iec if (grd%calving(i,j).ne.0.) grd%stored_heat(i,j)= & ! Need units of J sum(grd%stored_ice(i,j,:)) & ! initial stored ice in kg *grd%calving_hflx(i,j)*grd%area(i,j) & ! J/s/m2 x m^2 = J/s /grd%calving(i,j) ! /calving in kg/s enddo; enddo bergs%stored_heat_start=sum( grd%stored_heat(grd%isc:grd%iec,grd%jsc:grd%jec) ) call mpp_sum( bergs%stored_heat_start ) if (mpp_pe().eq.mpp_root_pe()) write(*,'(a,es13.6,a)') & 'diamonds, accumulate_calving: initial stored heat=',bergs%stored_heat_start,' J' endif remaining_dist=1. do k=1, nclasses grd%stored_ice(:,:,k)=grd%stored_ice(:,:,k)+bergs%dt*grd%calving(:,:)*bergs%distribution(k) remaining_dist=remaining_dist-bergs%distribution(k) enddo if (remaining_dist.lt.0.) then write(stderr(),*) 'diamonds, accumulate_calving: sum(distribution)>1!!!',remaining_dist call error_mesg('diamonds, accumulate_calving', 'calving is OVER distributed!', WARNING) endif net_calving_used=sum( grd%calving(grd%isc:grd%iec,grd%jsc:grd%jec) )*(1.-remaining_dist) bergs%net_calving_used=bergs%net_calving_used+net_calving_used*bergs%dt ! Remove the calving accounted for by accumulation grd%calving(:,:)=grd%calving(:,:)*remaining_dist ! Do the same for heat (no separate classes needed) grd%tmp(:,:)=bergs%dt*grd%calving_hflx(:,:)*grd%area(:,:)*(1.-remaining_dist) bergs%net_incoming_calving_heat_used=bergs%net_incoming_calving_heat_used+sum( grd%tmp(grd%isc:grd%iec,grd%jsc:grd%jec) ) grd%stored_heat(:,:)=grd%stored_heat(:,:)+grd%tmp(:,:) ! +=bergs%dt*grd%calving_hflx(:,:)*grd%area(:,:)*(1.-remaining_dist) grd%calving_hflx(:,:)=grd%calving_hflx(:,:)*remaining_dist end subroutine accumulate_calving ! ############################################################################## subroutine calve_icebergs(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables type(icebergs_gridded), pointer :: grd integer :: i,j,k,icnt,icntmax type(iceberg) :: newberg logical :: lret real :: xi, yj, ddt, calving_to_bergs, calved_to_berg, heat_to_bergs, heat_to_berg ! For convenience grd=>bergs%grd grd%real_calving(:,:,:)=0. calving_to_bergs=0. heat_to_bergs=0. icntmax=0 do k=1, nclasses do j=grd%jsc, grd%jec do i=grd%isc, grd%iec ddt=0.; icnt=0 do while (grd%stored_ice(i,j,k).ge.bergs%initial_mass(k)*bergs%mass_scaling(k)) newberg%lon=0.25*((grd%lon(i,j)+grd%lon(i-1,j-1))+(grd%lon(i-1,j)+grd%lon(i,j-1))) newberg%lat=0.25*((grd%lat(i,j)+grd%lat(i-1,j-1))+(grd%lat(i-1,j)+grd%lat(i,j-1))) !write(stderr(),*) 'diamonds, calve_icebergs: creating new iceberg at ',newberg%lon,newberg%lat lret=pos_within_cell(grd, newberg%lon, newberg%lat, i, j, xi, yj) if (.not.lret) then write(stderr(),*) 'diamonds, calve_icebergs: something went very wrong!',i,j,xi,yj call error_mesg('diamonds, calve_icebergs', 'berg is not in the correct cell!', FATAL) endif if (debug.and.(xi<0..or.xi>1..or.yj<0..or.yj>1.)) then write(stderr(),*) 'diamonds, calve_icebergs: something went very wrong!',i,j,xi,yj call error_mesg('diamonds, calve_icebergs', 'berg xi,yj is not correct!', FATAL) endif newberg%ine=i newberg%jne=j newberg%xi=xi newberg%yj=yj newberg%uvel=0. newberg%vvel=0. newberg%mass=bergs%initial_mass(k) newberg%thickness=bergs%initial_thickness(k) newberg%width=bergs%initial_width(k) newberg%length=bergs%initial_length(k) newberg%start_lon=newberg%lon newberg%start_lat=newberg%lat newberg%start_year=bergs%current_year newberg%start_day=bergs%current_yearday+ddt/86400. newberg%start_mass=bergs%initial_mass(k) newberg%mass_scaling=bergs%mass_scaling(k) newberg%mass_of_bits=0. newberg%heat_density=grd%stored_heat(i,j)/grd%stored_ice(i,j,k) ! This is in J/kg call add_new_berg_to_list(bergs%first, newberg) calved_to_berg=bergs%initial_mass(k)*bergs%mass_scaling(k) ! Units of kg ! Heat content heat_to_berg=calved_to_berg*newberg%heat_density ! Units of J grd%stored_heat(i,j)=grd%stored_heat(i,j)-heat_to_berg heat_to_bergs=heat_to_bergs+heat_to_berg ! Stored mass grd%stored_ice(i,j,k)=grd%stored_ice(i,j,k)-calved_to_berg calving_to_bergs=calving_to_bergs+calved_to_berg grd%real_calving(i,j,k)=grd%real_calving(i,j,k)+calved_to_berg/bergs%dt ddt=ddt+bergs%dt*2./17. ! Minor offset to start day icnt=icnt+1 bergs%nbergs_calved=bergs%nbergs_calved+1 bergs%nbergs_calved_by_class(k)=bergs%nbergs_calved_by_class(k)+1 enddo icntmax=max(icntmax,icnt) enddo enddo enddo if (debug.and.icntmax>1) write(stderr(),*) 'calve_icebergs: icnt=',icnt,' on',mpp_pe() bergs%net_calving_to_bergs=bergs%net_calving_to_bergs+calving_to_bergs bergs%net_heat_to_bergs=bergs%net_heat_to_bergs+heat_to_bergs end subroutine calve_icebergs ! ############################################################################## subroutine evolve_icebergs(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables type(icebergs_gridded), pointer :: grd real :: uvel1, vvel1, lon1, lat1, u1, v1, dxdl1, ax1, ay1 real :: uvel2, vvel2, lon2, lat2, u2, v2, dxdl2, ax2, ay2 real :: uvel3, vvel3, lon3, lat3, u3, v3, dxdl3, ax3, ay3 real :: uvel4, vvel4, lon4, lat4, u4, v4, dxdl4, ax4, ay4 real :: uveln, vveln, lonn, latn real :: x1, xdot1, xddot1, y1, ydot1, yddot1 real :: x2, xdot2, xddot2, y2, ydot2, yddot2 real :: x3, xdot3, xddot3, y3, ydot3, yddot3 real :: x4, xdot4, xddot4, y4, ydot4, yddot4 real :: xn, xdotn, yn, ydotn real :: r180_pi, dt, dt_2, dt_6, dydl, Rearth integer :: i, j integer :: i1,j1,i2,j2,i3,j3,i4,j4 real :: xi, yj logical :: bounced, on_tangential_plane, error_flag type(iceberg), pointer :: berg ! 4th order Runge-Kutta to solve: ! d/dt X = V, d/dt V = A ! with I.C.'s: ! X=X1 and V=V1 ! ! A1 = A(X1) ! X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1; A2=A(X2) ! X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3) ! X4 = X1+ dt*V3 ; V4 = V1+ dt*A3; A4=A(X4) ! ! Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 ! Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 ! For convenience grd=>bergs%grd ! Common constants r180_pi=1./pi_180 dt=bergs%dt dt_2=0.5*dt dt_6=dt/6. Rearth=6360.e3 berg=>bergs%first do while (associated(berg)) ! loop over all bergs if (.not. is_point_in_cell(bergs%grd, berg%lon, berg%lat, berg%ine, berg%jne) ) then write(stderr(),'(i4,a4,32i7)') mpp_pe(),'Lon',(i,i=grd%isd,grd%ied) do j=grd%jed,grd%jsd,-1 write(stderr(),'(2i4,32f7.1)') mpp_pe(),j,(grd%lon(i,j),i=grd%isd,grd%ied) enddo write(stderr(),'(i4,a4,32i7)') mpp_pe(),'Lat',(i,i=grd%isd,grd%ied) do j=grd%jed,grd%jsd,-1 write(stderr(),'(2i4,32f7.1)') mpp_pe(),j,(grd%lat(i,j),i=grd%isd,grd%ied) enddo call print_berg(stderr(), berg, 'evolve_iceberg, berg is not in proper starting cell') write(stderr(),'(a,i3,2(i4,3f8.2))') 'evolve_iceberg: pe,lon/lat(i,j)=', mpp_pe(), & berg%ine,berg%lon,grd%lon(berg%ine-1,berg%jne-1),grd%lon(berg%ine,berg%jne), & berg%jne,berg%lat,grd%lat(berg%ine-1,berg%jne-1),grd%lat(berg%ine,berg%jne) if (debug) call error_mesg('diamonds, evolve_iceberg','berg is in wrong starting cell!',FATAL) endif if (debug) call check_position(grd, berg, 'evolve_iceberg (top)') i=berg%ine j=berg%jne xi=berg%xi yj=berg%yj bounced=.false. on_tangential_plane=.false. if (berg%lat>89.) on_tangential_plane=.true. i1=i;j1=j if (bergs%add_weight_to_ocean .and. bergs%time_average_weight) & call spread_mass_across_ocean_cells(grd, i, j, xi, yj, berg%mass, berg%mass_of_bits, 0.25*berg%mass_scaling) ! A1 = A(X1) lon1=berg%lon; lat1=berg%lat if (on_tangential_plane) call rotpos_to_tang(lon1,lat1,x1,y1) dxdl1=r180_pi/(Rearth*cos(lat1*pi_180)) dydl=r180_pi/Rearth uvel1=berg%uvel; vvel1=berg%vvel if (on_tangential_plane) call rotvec_to_tang(lon1,uvel1,vvel1,xdot1,ydot1) u1=uvel1*dxdl1; v1=vvel1*dydl call accel(bergs, berg, i, j, xi, yj, lat1, uvel1, vvel1, uvel1, vvel1, dt_2, ax1, ay1) if (on_tangential_plane) call rotvec_to_tang(lon1,ax1,ay1,xddot1,yddot1) ! X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1; A2=A(X2) !if (debug) write(stderr(),*) 'diamonds, evolve: x2=...' if (on_tangential_plane) then x2=x1+dt_2*xdot1; y2=y1+dt_2*ydot1 xdot2=xdot1+dt_2*xddot1; ydot2=ydot1+dt_2*yddot1 call rotpos_from_tang(x2,y2,lon2,lat2) call rotvec_from_tang(lon2,xdot2,ydot2,uvel2,vvel2) else lon2=lon1+dt_2*u1; lat2=lat1+dt_2*v1 uvel2=uvel1+dt_2*ax1; vvel2=vvel1+dt_2*ay1 endif i=i1;j=j1;xi=berg%xi;yj=berg%yj call adjust_index_and_ground(grd, lon2, lat2, uvel2, vvel2, i, j, xi, yj, bounced, error_flag) i2=i; j2=j if (bergs%add_weight_to_ocean .and. bergs%time_average_weight) & call spread_mass_across_ocean_cells(grd, i, j, xi, yj, berg%mass, berg%mass_of_bits, 0.25*berg%mass_scaling) ! if (bounced.and.on_tangential_plane) call rotpos_to_tang(lon2,lat2,x2,y2) if (.not.error_flag) then if (debug .and. .not. is_point_in_cell(bergs%grd, lon2, lat2, i, j)) error_flag=.true. endif if (error_flag) then call print_fld(grd, grd%msk, 'msk') call print_fld(grd, grd%ssh, 'ssh') call print_fld(grd, grd%sst, 'sst') call print_fld(grd, grd%hi, 'hi') write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: isd,isc,iec,ied=',grd%isd,grd%isc,grd%iec,grd%ied write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: jsd,jsc,jec,jed=',grd%jsd,grd%jsc,grd%jec,grd%jed write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: i1,i2,i=',i1,i2,i write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: j1,j2,j=',j1,j2,j write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: lon1,lon2=',lon1,lon2,berg%lon write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: lat1,lat2=',lat1,lat2,berg%lat write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: u1,u2,u0=',uvel1,uvel2,berg%uvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: v1,v2,v0=',vvel1,vvel2,berg%vvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* ax1,ax2=',dt*ax1,dt*ax2 write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* ay1,ay2=',dt*ay1,dt*ay2 write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* u1,u2,u0=',dt*uvel1,dt*uvel2,dt*berg%uvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* v1,v2,v0=',dt*vvel1,dt*vvel2,dt*berg%vvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* u1,u2 (deg)=',dt*u1,dt*u2 write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* v1,v2 (deg)=',dt*v1,dt*v2 write(stderr(),*) 'Acceleration terms for position 1' error_flag=pos_within_cell(grd, lon1, lat1, i1, j1, xi, yj) call accel(bergs, berg, i1, j1, xi, yj, lat1, uvel1, vvel1, uvel1, vvel1, dt_2, ax1, ay1, debug_flag=.true.) call print_berg(stderr(), berg, 'evolve_iceberg, out of position at 2') write(stderr(),'(a,i3,a,2i4,4f8.3)') 'pe=',mpp_pe(),'pos2 i,j,lon,lat,xi,yj=',i,j,lon2,lat2,xi,yj write(stderr(),'(a,i3,a,4f8.3)') 'pe=',mpp_pe(),'pos2 box=',grd%lon(i-1,j-1),grd%lon(i,j),grd%lat(i-1,j-1),grd%lat(i,j) bounced=is_point_in_cell(bergs%grd, lon2, lat2, i, j, explain=.true.) call error_mesg('diamonds, evolve_iceberg','berg is out of posn at 2!',FATAL) endif dxdl2=r180_pi/(Rearth*cos(lat2*pi_180)) u2=uvel2*dxdl2; v2=vvel2*dydl call accel(bergs, berg, i, j, xi, yj, lat2, uvel2, vvel2, uvel1, vvel1, dt_2, ax2, ay2) if (on_tangential_plane) call rotvec_to_tang(lon2,ax2,ay2,xddot2,yddot2) ! X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3) !if (debug) write(stderr(),*) 'diamonds, evolve: x3=...' if (on_tangential_plane) then x3=x1+dt_2*xdot2; y3=y1+dt_2*ydot2 xdot3=xdot1+dt_2*xddot2; ydot3=ydot1+dt_2*yddot2 call rotpos_from_tang(x3,y3,lon3,lat3) call rotvec_from_tang(lon3,xdot3,ydot3,uvel3,vvel3) else lon3=lon1+dt_2*u2; lat3=lat1+dt_2*v2 uvel3=uvel1+dt_2*ax2; vvel3=vvel1+dt_2*ay2 endif i=i1;j=j1;xi=berg%xi;yj=berg%yj call adjust_index_and_ground(grd, lon3, lat3, uvel3, vvel3, i, j, xi, yj, bounced, error_flag) i3=i; j3=j if (bergs%add_weight_to_ocean .and. bergs%time_average_weight) & call spread_mass_across_ocean_cells(grd, i, j, xi, yj, berg%mass, berg%mass_of_bits, 0.25*berg%mass_scaling) ! if (bounced.and.on_tangential_plane) call rotpos_to_tang(lon3,lat3,x3,y3) if (.not.error_flag) then if (debug .and. .not. is_point_in_cell(bergs%grd, lon3, lat3, i, j)) error_flag=.true. endif if (error_flag) then call print_fld(grd, grd%msk, 'msk') call print_fld(grd, grd%ssh, 'ssh') call print_fld(grd, grd%sst, 'sst') call print_fld(grd, grd%hi, 'hi') write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: isd,isc,iec,ied=',grd%isd,grd%isc,grd%iec,grd%ied write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: jsd,jsc,jec,jed=',grd%jsd,grd%jsc,grd%jec,grd%jed write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: i1,i2,i3,i=',i1,i2,i3,i write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: j1,j2,j3,j=',j1,j2,j3,j write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: lon1,lon2,lon3=',lon1,lon2,lon3,berg%lon write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: lat1,lat2,lat3=',lat1,lat2,lat3,berg%lat write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: u1,u2,u3,u0=',uvel1,uvel2,uvel3,berg%uvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: v1,v2,v3,v0=',vvel1,vvel2,vvel3,berg%vvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* ax1,ax2,ax3=',dt*ax1,dt*ax2,dt*ax3 write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* ay1,ay2,ay3=',dt*ay1,dt*ay2,dt*ay3 write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* u1,u2,u3,u0=',dt*uvel1,dt*uvel2,dt*uvel3,dt*berg%uvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* v1,v2,v3,v0=',dt*vvel1,dt*vvel2,dt*vvel3,dt*berg%vvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* u1,u2,u3 (deg)=',dt*u1,dt*u2,dt*u3 write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* v1,v2,v3 (deg)=',dt*v1,dt*v2,dt*v3 write(stderr(),*) 'Acceleration terms for position 1' error_flag=pos_within_cell(grd, lon1, lat1, i1, j1, xi, yj) call accel(bergs, berg, i1, j1, xi, yj, lat1, uvel1, vvel1, uvel1, vvel1, dt_2, ax1, ay1, debug_flag=.true.) write(stderr(),*) 'Acceleration terms for position 2' error_flag=pos_within_cell(grd, lon2, lat2, i2, j2, xi, yj) call accel(bergs, berg, i2, j2, xi, yj, lat2, uvel2, vvel2, uvel1, vvel1, dt_2, ax2, ay2, debug_flag=.true.) call print_berg(stderr(), berg, 'evolve_iceberg, out of position at 3') write(stderr(),'(a,i3,a,2i4,4f8.3)') 'pe=',mpp_pe(),'pos3 i,j,lon,lat,xi,yj=',i,j,lon3,lat3,xi,yj write(stderr(),'(a,i3,a,4f8.3)') 'pe=',mpp_pe(),'pos3 box=',grd%lon(i-1,j-1),grd%lon(i,j),grd%lat(i-1,j-1),grd%lat(i,j) bounced=is_point_in_cell(bergs%grd, lon2, lat2, i, j, explain=.true.) call error_mesg('diamonds, evolve_iceberg','berg is out of posn at 3!',FATAL) endif dxdl3=r180_pi/(Rearth*cos(lat3*pi_180)) u3=uvel3*dxdl3; v3=vvel3*dydl call accel(bergs, berg, i, j, xi, yj, lat3, uvel3, vvel3, uvel1, vvel1, dt, ax3, ay3) if (on_tangential_plane) call rotvec_to_tang(lon3,ax3,ay3,xddot3,yddot3) ! X4 = X1+dt*V3 ; V4 = V1+dt*A3; A4=A(X4) !if (debug) write(stderr(),*) 'diamonds, evolve: x4=...' if (on_tangential_plane) then x4=x1+dt*xdot3; y4=y1+dt*ydot3 xdot4=xdot1+dt*xddot3; ydot4=ydot1+dt*yddot3 call rotpos_from_tang(x4,y4,lon4,lat4) call rotvec_from_tang(lon4,xdot4,ydot4,uvel4,vvel4) else lon4=lon1+dt*u3; lat4=lat1+dt*v3 uvel4=uvel1+dt*ax3; vvel4=vvel1+dt*ay3 endif i=i1;j=j1;xi=berg%xi;yj=berg%yj call adjust_index_and_ground(grd, lon4, lat4, uvel4, vvel4, i, j, xi, yj, bounced, error_flag) i4=i; j4=j ! if (bounced.and.on_tangential_plane) call rotpos_to_tang(lon4,lat4,x4,y4) if (.not.error_flag) then if (debug .and. .not. is_point_in_cell(bergs%grd, lon4, lat4, i, j)) error_flag=.true. endif if (error_flag) then call print_fld(grd, grd%msk, 'msk') call print_fld(grd, grd%ssh, 'ssh') call print_fld(grd, grd%sst, 'sst') call print_fld(grd, grd%hi, 'hi') write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: isd,isc,iec,ied=',grd%isd,grd%isc,grd%iec,grd%ied write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: jsd,jsc,jec,jed=',grd%jsd,grd%jsc,grd%jec,grd%jed write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: i1,i2,i3,i4,i=',i1,i2,i3,i4,i write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: j1,j2,j3,j4,j=',j1,j2,j3,j4,j write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: lon1,lon2,lon3,lon4=',lon1,lon2,lon3,lon4,berg%lon write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: lat1,lat2,lat3,lat4=',lat1,lat2,lat3,lat4,berg%lat write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: u1,u2,u3,u4,u0=',uvel1,uvel2,uvel3,uvel4,berg%uvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: v1,v2,v3,v4,v0=',vvel1,vvel2,vvel3,vvel4,berg%vvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* ax1,ax2,ax3,ax4=',dt*ax1,dt*ax2,dt*ax3,dt*ax4 write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* ay1,ay2,ay3,ay4=',dt*ay1,dt*ay2,dt*ay3,dt*ay4 write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* u1,u2,u3,u4,u0=',dt*uvel1,dt*uvel2,dt*uvel3,dt*uvel4,dt*berg%uvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* v1,v2,v3,v4,v0=',dt*vvel1,dt*vvel2,dt*vvel3,dt*vvel4,dt*berg%vvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* u1,u2,u3,u4 (deg)=',dt*u1,dt*u2,dt*u3,dt*u4 write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* v1,v2,v3,v4 (deg)=',dt*v1,dt*v2,dt*v3,dt*v4 write(stderr(),*) 'Acceleration terms for position 1' error_flag=pos_within_cell(grd, lon1, lat1, i1, j1, xi, yj) call accel(bergs, berg, i1, j1, xi, yj, lat1, uvel1, vvel1, uvel1, vvel1, dt_2, ax1, ay1, debug_flag=.true.) write(stderr(),*) 'Acceleration terms for position 2' error_flag=pos_within_cell(grd, lon2, lat2, i2, j2, xi, yj) call accel(bergs, berg, i2, j2, xi, yj, lat2, uvel2, vvel2, uvel1, vvel1, dt_2, ax2, ay2, debug_flag=.true.) write(stderr(),*) 'Acceleration terms for position 3' error_flag=pos_within_cell(grd, lon3, lat3, i3, j3, xi, yj) call accel(bergs, berg, i3, j3, xi, yj, lat3, uvel3, vvel3, uvel1, vvel1, dt, ax3, ay3, debug_flag=.true.) call print_berg(stderr(), berg, 'evolve_iceberg, out of position at 4') write(stderr(),'(a,i3,a,2i4,4f8.3)') 'pe=',mpp_pe(),'pos4 i,j,lon,lat,xi,yj=',i,j,lon4,lat4,xi,yj write(stderr(),'(a,i3,a,4f8.3)') 'pe=',mpp_pe(),'pos4 box=',grd%lon(i-1,j-1),grd%lon(i,j),grd%lat(i-1,j-1),grd%lat(i,j) bounced=is_point_in_cell(bergs%grd, lon2, lat2, i, j, explain=.true.) call error_mesg('diamonds, evolve_iceberg','berg is out of posn at 4!',FATAL) endif dxdl4=r180_pi/(Rearth*cos(lat4*pi_180)) u4=uvel4*dxdl4; v4=vvel4*dydl call accel(bergs, berg, i, j, xi, yj, lat4, uvel4, vvel4, uvel1, vvel1, dt, ax4, ay4) if (on_tangential_plane) call rotvec_to_tang(lon4,ax4,ay4,xddot4,yddot4) ! Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 ! Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 if (on_tangential_plane) then xn=x1+dt_6*( (xdot1+xdot4)+2.*(xdot2+xdot3) ) yn=y1+dt_6*( (ydot1+ydot4)+2.*(ydot2+ydot3) ) xdotn=xdot1+dt_6*( (xddot1+xddot4)+2.*(xddot2+xddot3) ) ydotn=ydot1+dt_6*( (yddot1+yddot4)+2.*(yddot2+yddot3) ) call rotpos_from_tang(xn,yn,lonn,latn) call rotvec_from_tang(lonn,xdotn,ydotn,uveln,vveln) else lonn=berg%lon+dt_6*( (u1+u4)+2.*(u2+u3) ) latn=berg%lat+dt_6*( (v1+v4)+2.*(v2+v3) ) uveln=berg%uvel+dt_6*( (ax1+ax4)+2.*(ax2+ax3) ) vveln=berg%vvel+dt_6*( (ay1+ay4)+2.*(ay2+ay3) ) endif i=i1;j=j1;xi=berg%xi;yj=berg%yj call adjust_index_and_ground(grd, lonn, latn, uveln, vveln, i, j, xi, yj, bounced, error_flag) if (bergs%add_weight_to_ocean .and. bergs%time_average_weight) & call spread_mass_across_ocean_cells(grd, i, j, xi, yj, berg%mass, berg%mass_of_bits, 0.25*berg%mass_scaling) if (.not.error_flag) then if (.not. is_point_in_cell(bergs%grd, lonn, latn, i, j)) error_flag=.true. endif if (error_flag) then call print_fld(grd, grd%msk, 'msk') call print_fld(grd, grd%ssh, 'ssh') call print_fld(grd, grd%sst, 'sst') call print_fld(grd, grd%hi, 'hi') write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: isd,isc,iec,ied=',grd%isd,grd%isc,grd%iec,grd%ied write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: jsd,jsc,jec,jed=',grd%jsd,grd%jsc,grd%jec,grd%jed write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: i1,i2,i3,i4,i=',i1,i2,i3,i4,i write(stderr(),'(a,6i5)') 'diamonds, evolve_iceberg: j1,j2,j3,j4,j=',j1,j2,j3,j4,j write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: lon1,lon2,lon3,lon4,lonn=',lon1,lon2,lon3,lon4,lonn,berg%lon write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: lat1,lat2,lat3,lat4,latn=',lat1,lat2,lat3,lat4,latn,berg%lat write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: u1,u2,u3,u4,un,u0=',uvel1,uvel2,uvel3,uvel4,uveln,berg%uvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: v1,v2,v3,v4,vn,v0=',vvel1,vvel2,vvel3,vvel4,vveln,berg%vvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* ax1,ax2,ax3,ax4,axn=',dt*ax1,dt*ax2,dt*ax3,dt*ax4,dt_6*( (ax1+ax4)+2.*(ax2+ax3) ) write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* ay1,ay2,ay3,ay4,ayn=',dt*ay1,dt*ay2,dt*ay3,dt*ay4,dt_6*( (ay1+ay4)+2.*(ay2+ay3) ) write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* u1,u2,u3,u4,un,u0=',dt*uvel1,dt*uvel2,dt*uvel3,dt*uvel4,dt*uveln,dt*berg%uvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* v1,v2,v3,v4,vn,v0=',dt*vvel1,dt*vvel2,dt*vvel3,dt*vvel4,dt*vveln,dt*berg%vvel write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* u1,u2,u3,u4,u_rk (deg)=',dt*u1,dt*u2,dt*u3,dt*u4,dt_6*( (u1+u4)+2.*(u2+u3) ) write(stderr(),'(a,6es)') 'diamonds, evolve_iceberg: dt* v1,v2,v3,v4,v_rk (deg)=',dt*v1,dt*v2,dt*v3,dt*v4,dt_6*( (v1+v4)+2.*(v2+v3) ) write(stderr(),*) 'diamonds, evolve_iceberg: on_tangential_plane=',on_tangential_plane write(stderr(),*) 'Acceleration terms for position 1' error_flag=pos_within_cell(grd, lon1, lat1, i1, j1, xi, yj) call accel(bergs, berg, i1, j1, xi, yj, lat1, uvel1, vvel1, uvel1, vvel1, dt_2, ax1, ay1, debug_flag=.true.) write(stderr(),*) 'Acceleration terms for position 2' error_flag=pos_within_cell(grd, lon2, lat2, i2, j2, xi, yj) call accel(bergs, berg, i2, j2, xi, yj, lat2, uvel2, vvel2, uvel1, vvel1, dt_2, ax2, ay2, debug_flag=.true.) write(stderr(),*) 'Acceleration terms for position 3' error_flag=pos_within_cell(grd, lon3, lat3, i3, j3, xi, yj) call accel(bergs, berg, i3, j3, xi, yj, lat3, uvel3, vvel3, uvel1, vvel1, dt, ax3, ay3, debug_flag=.true.) write(stderr(),*) 'Acceleration terms for position 4' error_flag=pos_within_cell(grd, lon4, lat4, i4, j4, xi, yj) call accel(bergs, berg, i4, j4, xi, yj, lat4, uvel4, vvel4, uvel1, vvel1, dt, ax4, ay4, debug_flag=.true.) write(stderr(),'(a,i3,a,2i4,4f8.3)') 'pe=',mpp_pe(),'posn i,j,lon,lat,xi,yj=',i,j,lonn,latn,xi,yj write(stderr(),'(a,i3,a,4f8.3)') 'pe=',mpp_pe(),'posn box=',grd%lon(i-1,j-1),grd%lon(i,j),grd%lat(i-1,j-1),grd%lat(i,j) call print_berg(stderr(), berg, 'evolve_iceberg, out of cell at end!') bounced=is_point_in_cell(bergs%grd, lonn, latn, i, j, explain=.true.) if (debug) call error_mesg('diamonds, evolve_iceberg','berg is out of posn at end!',FATAL) write(stderr(),'(i4,a4,32i7)') mpp_pe(),'Lon',(i,i=grd%isd,grd%ied) do j=grd%jed,grd%jsd,-1 write(stderr(),'(2i4,32f7.1)') mpp_pe(),j,(grd%lon(i,j),i=grd%isd,grd%ied) enddo write(stderr(),'(i4,a4,32i7)') mpp_pe(),'Lat',(i,i=grd%isd,grd%ied) do j=grd%jed,grd%jsd,-1 write(stderr(),'(2i4,32f7.1)') mpp_pe(),j,(grd%lat(i,j),i=grd%isd,grd%ied) enddo endif berg%lon=lonn berg%lat=latn berg%uvel=uveln berg%vvel=vveln berg%ine=i berg%jne=j berg%xi=xi berg%yj=yj !call interp_flds(grd, i, j, xi, yj, berg%uo, berg%vo, berg%ui, berg%vi, berg%ua, berg%va, berg%ssh_x, berg%ssh_y, berg%sst) !if (debug) call print_berg(stderr(), berg, 'evolve_iceberg, final posn.') if (debug) call check_position(grd, berg, 'evolve_iceberg (bot)') berg=>berg%next enddo ! loop over all bergs contains subroutine rotpos_to_tang(lon, lat, x, y) ! Arguments real, intent(in) :: lon, lat real, intent(out) :: x, y ! Local variables real :: r,colat,clon,slon if (lat>90.) then write(stderr(),*) 'diamonds, rotpos_to_tang: lat>90 already!',lat call error_mesg('diamonds, rotpos_to_tang','Something went very wrong!',FATAL) endif if (lat==90.) then write(stderr(),*) 'diamonds, rotpos_to_tang: lat==90 already!',lat call error_mesg('diamonds, rotpos_to_tang','Something went wrong!',FATAL) endif colat=90.-lat r=Rearth*(colat*pi_180) clon=cos(lon*pi_180) slon=sin(lon*pi_180) x=r*clon y=r*slon end subroutine rotpos_to_tang subroutine rotpos_from_tang(x, y, lon, lat) ! Arguments real, intent(in) :: x, y real, intent(out) :: lon, lat ! Local variables real :: r r=sqrt(x**2+y**2) lat=90.-(r180_pi*r/Rearth) lon=r180_pi*acos(x/r)*sign(1.,y) end subroutine rotpos_from_tang subroutine rotvec_to_tang(lon, uvel, vvel, xdot, ydot) ! Arguments real, intent(in) :: lon, uvel, vvel real, intent(out) :: xdot, ydot ! Local variables real :: clon,slon clon=cos(lon*pi_180) slon=sin(lon*pi_180) xdot=-slon*uvel-clon*vvel ydot=clon*uvel-slon*vvel end subroutine rotvec_to_tang subroutine rotvec_from_tang(lon, xdot, ydot, uvel, vvel) ! Arguments real, intent(in) :: lon, xdot, ydot real, intent(out) :: uvel, vvel ! Local variables real :: clon,slon clon=cos(lon*pi_180) slon=sin(lon*pi_180) uvel=-slon*xdot+clon*ydot vvel=-clon*xdot-slon*ydot end subroutine rotvec_from_tang ! ############################################################################## subroutine adjust_index_and_ground(grd, lon, lat, uvel, vvel, i, j, xi, yj, bounced, error) ! Arguments type(icebergs_gridded), pointer :: grd real, intent(inout) :: lon, lat, uvel, vvel, xi, yj integer, intent(inout) :: i,j logical, intent(out) :: bounced, error ! Local variables logical lret, lpos real, parameter :: posn_eps=0.05 integer :: icount, i0, j0, inm, jnm real :: xi0, yj0, lon0, lat0 bounced=.false. error=.false. lon0=lon; lat0=lat ! original position i0=i; j0=j ! original i,j lret=pos_within_cell(grd, lon, lat, i, j, xi, yj) xi0=xi; yj0=yj ! original xi,yj if (debug) then !Sanity check lret, xi and yj lret=is_point_in_cell(grd, lon, lat, i, j) if (xi<0. .or. xi>1. .or. yj<0. .or. yj>1.) then if (lret) then write(stderr(),*) 'diamonds, adjust: WARNING!!! lret=T but |xi,yj|>1',mpp_pe() write(stderr(),*) 'diamonds, adjust: xi=',xi,' lon=',lon write(stderr(),*) 'diamonds, adjust: x3 x2=',grd%lon(i-1,j),grd%lon(i,j) write(stderr(),*) 'diamonds, adjust: x0 x1=',grd%lon(i-1,j-1),grd%lon(i,j-1) write(stderr(),*) 'diamonds, adjust: yi=',yj,' lat=',lat write(stderr(),*) 'diamonds, adjust: y3 y2=',grd%lat(i-1,j),grd%lat(i,j) write(stderr(),*) 'diamonds, adjust: y0 y1=',grd%lat(i-1,j-1),grd%lat(i,j-1) lret=is_point_in_cell(grd, lon, lat, i, j, explain=.true.) write(stderr(),*) 'diamonds, adjust: fn is_point_in_cell=',lret lret=pos_within_cell(grd, lon, lat, i, j, xi, yj, explain=.true.) write(stderr(),*) 'diamonds, adjust: fn pos_within_cell=',lret write(0,*) 'This should never happen!' error=.true.; return endif else if (.not.lret) then write(stderr(),*) 'diamonds, adjust: WARNING!!! lret=F but |xi,yj|<1',mpp_pe() write(stderr(),*) 'diamonds, adjust: xi=',xi,' lon=',lon write(stderr(),*) 'diamonds, adjust: x3 x2=',grd%lon(i-1,j),grd%lon(i,j) write(stderr(),*) 'diamonds, adjust: x0 x1=',grd%lon(i-1,j-1),grd%lon(i,j-1) write(stderr(),*) 'diamonds, adjust: yi=',yj,' lat=',lat write(stderr(),*) 'diamonds, adjust: y3 y2=',grd%lat(i-1,j),grd%lat(i,j) write(stderr(),*) 'diamonds, adjust: y0 y1=',grd%lat(i-1,j-1),grd%lat(i,j-1) lret=is_point_in_cell(grd, lon, lat, i, j, explain=.true.) write(stderr(),*) 'diamonds, adjust: fn is_point_in_cell=',lret lret=pos_within_cell(grd, lon, lat, i, j, xi, yj, explain=.true.) write(stderr(),*) 'diamonds, adjust: fn pos_within_cell=',lret write(0,*) 'This should never happen!' error=.true.; return endif endif lret=pos_within_cell(grd, lon, lat, i, j, xi, yj) endif ! debug if (lret) return ! Berg was already in cell icount=0 inm=i0; jnm=j0 ! original i,j do while ( .not.lret.and. icount<4 ) icount=icount+1 if (xi.lt.0.) then if (inm>grd%isd) then inm=inm-1 endif elseif (xi.gt.1.) then if (inmgrd%jsd) then jnm=jnm-1 endif elseif (yj.gt.1.) then if (jnm1) then write(stderr(),*) 'pe=',mpp_pe(),'diamonds, adjust: inm,i0,inm-i0=',inm,i0,inm-i0 !stop 'Moved too far in i without mask!' endif if (abs(jnm-j0)>1) then write(stderr(),*) 'pe=',mpp_pe(),'diamonds, adjust: jnm,i0,jnm-j0=',jnm,j0,inm-j0 !stop 'Moved too far in j without mask!' endif lret=pos_within_cell(grd, lon, lat, i0, j0, xi, yj) icount=0 do while ( .not.lret.and. icount<4 ) icount=icount+1 if (xi.lt.0.) then if (i>grd%isd) then if (grd%msk(i-1,j)>0.) then if (i>grd%isd+1) i=i-1 else !write(stderr(),'(a,6f8.3,i)') 'diamonds, adjust: bouncing berg from west',lon,lat,xi,yj,uvel,vvel,mpp_pe() bounced=.true. endif endif elseif (xi.gt.1.) then if (i0.) then if (igrd%jsd) then if (grd%msk(i,j-1)>0.) then if (j>grd%jsd+1) j=j-1 else !write(stderr(),'(a,6f8.3,i)') 'diamonds, adjust: bouncing berg from south',lon,lat,xi,yj,uvel,vvel,mpp_pe() bounced=.true. endif endif elseif (yj.gt.1.) then if (j0.) then if (j1.) xi=1.-posn_eps if (xi<0.) xi=posn_eps if (yj>1.) yj=1.-posn_eps if (yj<0.) yj=posn_eps lon=bilin(grd, grd%lon, i, j, xi, yj) lat=bilin(grd, grd%lat, i, j, xi, yj) endif if (debug) then if (grd%msk(i,j)==0.) stop 'diamonds, adjust: Berg is in land! This should not happen...' endif lret=pos_within_cell(grd, lon, lat, i, j, xi, yj) ! Update xi and yj enddo !if (debug) then ! if (abs(i-i0)>2) then ! stop 'diamonds, adjust: Moved too far in i!' ! endif ! if (abs(j-j0)>2) then ! stop 'diamonds, adjust: Moved too far in j!' ! endif !endif if (.not.bounced.and.lret.and.grd%msk(i,j)>0.) return ! Landed in ocean without bouncing so all is well if (.not.bounced.and..not.lret) then ! This implies the berg traveled many cells if (debug) then write(stderr(),*) 'diamonds, adjust: lon0, lat0=',lon0,lat0 write(stderr(),*) 'diamonds, adjust: xi0, yj0=',xi0,yj0 write(stderr(),*) 'diamonds, adjust: i0,j0=',i0,j0 write(stderr(),*) 'diamonds, adjust: lon, lat=',lon,lat write(stderr(),*) 'diamonds, adjust: xi,yj=',xi,yj write(stderr(),*) 'diamonds, adjust: i,j=',i,j write(stderr(),*) 'diamonds, adjust: inm,jnm=',inm,jnm write(stderr(),*) 'diamonds, adjust: icount=',icount lret=pos_within_cell(grd, lon, lat, i, j, xi, yj, explain=.true.) write(stderr(),*) 'diamonds, adjust: lret=',lret endif call error_mesg('diamonds, adjust', 'Berg iterated many times without bouncing!', WARNING) endif if (xi>1.) xi=1.-posn_eps if (xi<0.) xi=posn_eps if (yj>1.) yj=1.-posn_eps if (yj<0.) yj=posn_eps lon=bilin(grd, grd%lon, i, j, xi, yj) lat=bilin(grd, grd%lat, i, j, xi, yj) lret=pos_within_cell(grd, lon, lat, i, j, xi, yj) ! Update xi and yj if (.not. lret) then write(stderr(),*) 'diamonds, adjust: Should not get here! Berg is not in cell after adjustment' if (debug) error=.true. endif end subroutine adjust_index_and_ground end subroutine evolve_icebergs ! ############################################################################## subroutine send_bergs_to_other_pes(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables type(iceberg), pointer :: kick_the_bucket, this integer :: nbergs_to_send_e, nbergs_to_send_w integer :: nbergs_to_send_n, nbergs_to_send_s integer :: nbergs_rcvd_from_e, nbergs_rcvd_from_w integer :: nbergs_rcvd_from_n, nbergs_rcvd_from_s integer, parameter :: buffer_width=18 type(icebergs_gridded), pointer :: grd integer :: i, nbergs_start, nbergs_end ! For convenience grd=>bergs%grd if (debug) then nbergs_start=count_bergs(bergs) endif ! Find number of bergs that headed east/west nbergs_to_send_e=0 nbergs_to_send_w=0 if (associated(bergs%first)) then this=>bergs%first do while (associated(this)) if (this%ine.gt.bergs%grd%iec) then kick_the_bucket=>this this=>this%next nbergs_to_send_e=nbergs_to_send_e+1 call pack_berg_into_buffer2(kick_the_bucket, bergs%obuffer_e, nbergs_to_send_e) call move_trajectory(bergs, kick_the_bucket) call delete_iceberg_from_list(bergs%first,kick_the_bucket) elseif (this%ine.lt.bergs%grd%isc) then kick_the_bucket=>this this=>this%next nbergs_to_send_w=nbergs_to_send_w+1 call pack_berg_into_buffer2(kick_the_bucket, bergs%obuffer_w, nbergs_to_send_w) call move_trajectory(bergs, kick_the_bucket) call delete_iceberg_from_list(bergs%first,kick_the_bucket) else this=>this%next endif enddo endif ! Send bergs east if (grd%pe_E.ne.NULL_PE) then call mpp_send(nbergs_to_send_e, plen=1, to_pe=grd%pe_E) if (nbergs_to_send_e.gt.0) then call mpp_send(bergs%obuffer_e%data, nbergs_to_send_e*buffer_width, grd%pe_E) endif endif ! Send bergs west if (grd%pe_W.ne.NULL_PE) then call mpp_send(nbergs_to_send_w, plen=1, to_pe=grd%pe_W) if (nbergs_to_send_w.gt.0) then call mpp_send(bergs%obuffer_w%data, nbergs_to_send_w*buffer_width, grd%pe_W) endif endif ! Receive bergs from west if (grd%pe_W.ne.NULL_PE) then nbergs_rcvd_from_w=-999 call mpp_recv(nbergs_rcvd_from_w, glen=1, from_pe=grd%pe_W) if (nbergs_rcvd_from_w.lt.0) then write(stderr(),*) 'pe=',mpp_pe(),' received a bad number',nbergs_rcvd_from_w,' from',grd%pe_W,' (W) !!!!!!!!!!!!!!!!!!!!!!' endif if (nbergs_rcvd_from_w.gt.0) then call increase_ibuffer(bergs%ibuffer_w, nbergs_rcvd_from_w) call mpp_recv(bergs%ibuffer_w%data, nbergs_rcvd_from_w*buffer_width, grd%pe_W) do i=1, nbergs_rcvd_from_w call unpack_berg_from_buffer2(bergs%first, bergs%ibuffer_w, i) enddo endif else nbergs_rcvd_from_w=0 endif ! Receive bergs from east if (grd%pe_E.ne.NULL_PE) then nbergs_rcvd_from_e=-999 call mpp_recv(nbergs_rcvd_from_e, glen=1, from_pe=grd%pe_E) if (nbergs_rcvd_from_e.lt.0) then write(stderr(),*) 'pe=',mpp_pe(),' received a bad number',nbergs_rcvd_from_e,' from',grd%pe_E,' (E) !!!!!!!!!!!!!!!!!!!!!!' endif if (nbergs_rcvd_from_e.gt.0) then call increase_ibuffer(bergs%ibuffer_e, nbergs_rcvd_from_e) call mpp_recv(bergs%ibuffer_e%data, nbergs_rcvd_from_e*buffer_width, grd%pe_E) do i=1, nbergs_rcvd_from_e call unpack_berg_from_buffer2(bergs%first, bergs%ibuffer_e, i) enddo endif else nbergs_rcvd_from_e=0 endif ! Find number of bergs that headed north/south ! (note: this block should technically go ahead of the E/W recv block above ! to handle arbitrary orientation of PEs. But for simplicity, it is ! here to accomodate diagonal transfer of bergs between PEs -AJA) nbergs_to_send_n=0 nbergs_to_send_s=0 if (associated(bergs%first)) then this=>bergs%first do while (associated(this)) if (this%jne.gt.bergs%grd%jec) then kick_the_bucket=>this this=>this%next nbergs_to_send_n=nbergs_to_send_n+1 call pack_berg_into_buffer2(kick_the_bucket, bergs%obuffer_n, nbergs_to_send_n) call move_trajectory(bergs, kick_the_bucket) call delete_iceberg_from_list(bergs%first,kick_the_bucket) elseif (this%jne.lt.bergs%grd%jsc) then kick_the_bucket=>this this=>this%next nbergs_to_send_s=nbergs_to_send_s+1 call pack_berg_into_buffer2(kick_the_bucket, bergs%obuffer_s, nbergs_to_send_s) call move_trajectory(bergs, kick_the_bucket) call delete_iceberg_from_list(bergs%first,kick_the_bucket) else this=>this%next endif enddo endif ! Send bergs north if (grd%pe_N.ne.NULL_PE) then call mpp_send(nbergs_to_send_n, plen=1, to_pe=grd%pe_N) if (nbergs_to_send_n.gt.0) then call mpp_send(bergs%obuffer_n%data, nbergs_to_send_n*buffer_width, grd%pe_N) endif endif ! Send bergs south if (grd%pe_S.ne.NULL_PE) then call mpp_send(nbergs_to_send_s, plen=1, to_pe=grd%pe_S) if (nbergs_to_send_s.gt.0) then call mpp_send(bergs%obuffer_s%data, nbergs_to_send_s*buffer_width, grd%pe_S) endif endif ! Receive bergs from south if (grd%pe_S.ne.NULL_PE) then nbergs_rcvd_from_s=-999 call mpp_recv(nbergs_rcvd_from_s, glen=1, from_pe=grd%pe_S) if (nbergs_rcvd_from_s.lt.0) then write(stderr(),*) 'pe=',mpp_pe(),' received a bad number',nbergs_rcvd_from_s,' from',grd%pe_S,' (S) !!!!!!!!!!!!!!!!!!!!!!' endif if (nbergs_rcvd_from_s.gt.0) then call increase_ibuffer(bergs%ibuffer_s, nbergs_rcvd_from_s) call mpp_recv(bergs%ibuffer_s%data, nbergs_rcvd_from_s*buffer_width, grd%pe_S) do i=1, nbergs_rcvd_from_s call unpack_berg_from_buffer2(bergs%first, bergs%ibuffer_s, i) enddo endif else nbergs_rcvd_from_s=0 endif ! Receive bergs from north if (grd%pe_N.ne.NULL_PE) then nbergs_rcvd_from_n=-999 call mpp_recv(nbergs_rcvd_from_n, glen=1, from_pe=grd%pe_N) if (nbergs_rcvd_from_n.lt.0) then write(stderr(),*) 'pe=',mpp_pe(),' received a bad number',nbergs_rcvd_from_n,' from',grd%pe_N,' (N) !!!!!!!!!!!!!!!!!!!!!!' endif if (nbergs_rcvd_from_n.gt.0) then call increase_ibuffer(bergs%ibuffer_n, nbergs_rcvd_from_n) call mpp_recv(bergs%ibuffer_n%data, nbergs_rcvd_from_n*buffer_width, grd%pe_N) do i=1, nbergs_rcvd_from_n call unpack_berg_from_buffer2(bergs%first, bergs%ibuffer_n, i) enddo endif else nbergs_rcvd_from_n=0 endif if (debug) then nbergs_end=count_bergs(bergs) i=nbergs_rcvd_from_n+nbergs_rcvd_from_s+nbergs_rcvd_from_e+nbergs_rcvd_from_w & -nbergs_to_send_n-nbergs_to_send_s-nbergs_to_send_e-nbergs_to_send_w if (nbergs_end-(nbergs_start+i).ne.0) then write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_end=',nbergs_end,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_start=',nbergs_start,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: delta=',i,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: error=',nbergs_end-(nbergs_start+i),' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_to_send_n=',nbergs_to_send_n,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_to_send_s=',nbergs_to_send_s,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_to_send_e=',nbergs_to_send_e,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_to_send_w=',nbergs_to_send_w,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_rcvd_from_n=',nbergs_rcvd_from_n,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_rcvd_from_s=',nbergs_rcvd_from_s,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_rcvd_from_e=',nbergs_rcvd_from_e,' on PE',mpp_pe() write(stderr(),'(a,i,a,i4)') 'diamonds, send_bergs_to_other_pes: nbergs_rcvd_from_w=',nbergs_rcvd_from_w,' on PE',mpp_pe() call error_mesg('diamonds, send_bergs_to_other_pes:', 'We lost some bergs!', FATAL) endif endif if (debug) then i=0 this=>bergs%first do while (associated(this)) call check_position(grd, this, 'exchange (bot)') if (this%ine.lt.bergs%grd%isc .or. & this%ine.gt.bergs%grd%iec .or. & this%jne.lt.bergs%grd%jsc .or. & this%jne.gt.bergs%grd%jec) i=i+1 this=>this%next enddo ! while call mpp_sum(i) if (i>0 .and. mpp_pe()==mpp_root_pe()) then write(stderr(),'(a,i4)') 'diamonds, send_bergs_to_other_pes: # of bergs outside computational domain = ',i call error_mesg('diamonds, send_bergs_to_other_pes:', 'there are bergs still in halos!', FATAL) endif ! root_pe endif ! debug call mpp_sync_self() contains subroutine pack_berg_into_buffer2(berg, buff, n) ! Arguments type(iceberg), pointer :: berg type(buffer), pointer :: buff integer, intent(in) :: n ! Local variables if (.not.associated(buff)) call increase_buffer(buff,delta_buf) if (n>buff%size) call increase_buffer(buff,delta_buf) buff%data(1,n)=berg%lon buff%data(2,n)=berg%lat buff%data(3,n)=berg%uvel buff%data(4,n)=berg%vvel buff%data(5,n)=berg%xi buff%data(6,n)=berg%yj buff%data(7,n)=berg%start_lon buff%data(8,n)=berg%start_lat buff%data(9,n)=float(berg%start_year) buff%data(10,n)=berg%start_day buff%data(11,n)=berg%start_mass buff%data(12,n)=berg%mass buff%data(13,n)=berg%thickness buff%data(14,n)=berg%width buff%data(15,n)=berg%length buff%data(16,n)=berg%mass_scaling buff%data(17,n)=berg%mass_of_bits buff%data(18,n)=berg%heat_density end subroutine pack_berg_into_buffer2 subroutine increase_buffer(old,delta) ! Arguments type(buffer), pointer :: old integer, intent(in) :: delta ! Local variables type(buffer), pointer :: new integer :: new_size if (.not.associated(old)) then new_size=delta else new_size=old%size+delta endif allocate(new) allocate(new%data(buffer_width,new_size)) new%size=new_size if (associated(old)) then new%data(:,1:old%size)=old%data(:,1:old%size) deallocate(old%data) deallocate(old) endif old=>new !write(stderr(),*) 'diamonds, increase_buffer',mpp_pe(),' increased to',new_size end subroutine increase_buffer subroutine unpack_berg_from_buffer2(first, buff, n) ! Arguments type(iceberg), pointer :: first type(buffer), pointer :: buff integer, intent(in) :: n ! Local variables !real :: lon, lat, uvel, vvel, xi, yj !real :: start_lon, start_lat, start_day, start_mass !integer :: ine, jne, start_year logical :: lres type(iceberg) :: localberg localberg%lon=buff%data(1,n) localberg%lat=buff%data(2,n) localberg%uvel=buff%data(3,n) localberg%vvel=buff%data(4,n) localberg%xi=buff%data(5,n) localberg%yj=buff%data(6,n) localberg%start_lon=buff%data(7,n) localberg%start_lat=buff%data(8,n) localberg%start_year=int(buff%data(9,n)+0.5) localberg%start_day=buff%data(10,n) localberg%start_mass=buff%data(11,n) localberg%mass=buff%data(12,n) localberg%thickness=buff%data(13,n) localberg%width=buff%data(14,n) localberg%length=buff%data(15,n) localberg%mass_scaling=buff%data(16,n) localberg%mass_of_bits=buff%data(17,n) localberg%heat_density=buff%data(18,n) lres=find_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne) if (lres) then lres=pos_within_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne, localberg%xi, localberg%yj) call add_new_berg_to_list(first, localberg) else lres=find_cell_wide(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne) if (lres) then lres=pos_within_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne, localberg%xi, localberg%yj) call add_new_berg_to_list(first, localberg) else write(stderr(),'("diamonds, unpack_berg_from_buffer pe=(",i3,a,2i4,a,2f8.2)') mpp_pe(),') Failed to find i,j=',localberg%ine,localberg%jne,' for lon,lat=',localberg%lon,localberg%lat write(stderr(),*) localberg%lon,localberg%lat write(stderr(),*) localberg%uvel,localberg%vvel write(stderr(),*) grd%isc,grd%iec,grd%jsc,grd%jec write(stderr(),*) grd%isd,grd%ied,grd%jsd,grd%jed write(stderr(),*) grd%lon(grd%isc-1,grd%jsc-1),grd%lon(grd%iec,grd%jsc) write(stderr(),*) grd%lat(grd%isc-1,grd%jsc-1),grd%lat(grd%iec,grd%jec) write(stderr(),*) grd%lon(grd%isd,grd%jsd),grd%lon(grd%ied,grd%jsd) write(stderr(),*) grd%lat(grd%isd,grd%jsd),grd%lat(grd%ied,grd%jed) write(stderr(),*) lres call error_mesg('diamonds, unpack_berg_from_buffer', 'can not find a cell to place berg in!', FATAL) endif endif end subroutine unpack_berg_from_buffer2 subroutine increase_ibuffer(old,delta) ! Arguments type(buffer), pointer :: old integer, intent(in) :: delta ! Local variables type(buffer), pointer :: new integer :: new_size, old_size if (.not.associated(old)) then new_size=delta+delta_buf old_size=0 else old_size=old%size if (deltanew !write(stderr(),*) 'diamonds, increase_ibuffer',mpp_pe(),' increased to',new_size endif end subroutine increase_ibuffer end subroutine send_bergs_to_other_pes ! ############################################################################## subroutine icebergs_init(bergs, & gni, gnj, layout, io_layout, axes, maskmap, x_cyclic, tripolar_grid, & dt, Time, ice_lon, ice_lat, ice_wet, ice_dx, ice_dy, ice_area, & cos_rot, sin_rot) ! Arguments type(icebergs), pointer :: bergs integer, intent(in) :: gni, gnj, layout(2), io_layout(2), axes(2) logical, intent(in), optional :: maskmap(:,:) logical, intent(in) :: x_cyclic, tripolar_grid real, intent(in) :: dt type (time_type), intent(in) :: Time ! current time real, dimension(:,:), intent(in) :: ice_lon, ice_lat, ice_wet real, dimension(:,:), intent(in) :: ice_dx, ice_dy, ice_area real, dimension(:,:), intent(in) :: cos_rot, sin_rot ! Namelist parameters (and defaults) integer :: halo=4 ! Width of halo region integer :: traj_sample_hrs=24 ! Period between sampling of position for trajectory storage integer :: verbose_hrs=24 ! Period between verbose messages real :: rho_bergs=850. ! Density of icebergs real :: LoW_ratio=1.5 ! Initial ratio L/W for newly calved icebergs real :: bergy_bit_erosion_fraction=0. ! Fraction of erosion melt flux to divert to bergy bits real :: sicn_shift=0. ! Shift of sea-ice concentration in erosion flux modulation (0bergs%grd ! For convenience to avoid bergs%grd%X !write(stderrunit,*) 'diamonds: allocating domain' allocate(grd%domain) ! Clocks bergs%clock=mpp_clock_id( 'Icebergs', flags=clock_flag_default, grain=CLOCK_COMPONENT ) bergs%clock_mom=mpp_clock_id( 'Icebergs-momentum', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT ) bergs%clock_the=mpp_clock_id( 'Icebergs-thermodyn', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT ) bergs%clock_int=mpp_clock_id( 'Icebergs-interface', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT ) bergs%clock_cal=mpp_clock_id( 'Icebergs-calving', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT ) bergs%clock_com=mpp_clock_id( 'Icebergs-communication', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT ) bergs%clock_ini=mpp_clock_id( 'Icebergs-initialization', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT ) bergs%clock_ior=mpp_clock_id( 'Icebergs-I/O read', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT ) bergs%clock_iow=mpp_clock_id( 'Icebergs-I/O write', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT ) bergs%clock_dia=mpp_clock_id( 'Icebergs-diagnostics', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT ) call mpp_clock_begin(bergs%clock) call mpp_clock_begin(bergs%clock_ini) ! Set up iceberg domain !write(stderrunit,*) 'diamonds: defining domain' if(tripolar_grid) then call mpp_define_domains( (/1,gni,1,gnj/), layout, grd%domain, & ! maskmap=maskmap, & xflags=CYCLIC_GLOBAL_DOMAIN, xhalo=halo, & yflags=FOLD_NORTH_EDGE, yhalo=halo, name='diamond') else if(x_cyclic) then call mpp_define_domains( (/1,gni,1,gnj/), layout, grd%domain, & ! maskmap=maskmap, & xflags=CYCLIC_GLOBAL_DOMAIN, & xhalo=halo, yhalo=halo, name='diamond') else call mpp_define_domains( (/1,gni,1,gnj/), layout, grd%domain, & ! maskmap=maskmap, & xhalo=halo, yhalo=halo, name='diamond') endif call mpp_define_io_domain(grd%domain, io_layout) !write(stderrunit,*) 'diamond: get compute domain' call mpp_get_compute_domain( grd%domain, grd%isc, grd%iec, grd%jsc, grd%jec ) call mpp_get_data_domain( grd%domain, grd%isd, grd%ied, grd%jsd, grd%jed ) call mpp_get_neighbor_pe(grd%domain, NORTH, grd%pe_N) call mpp_get_neighbor_pe(grd%domain, SOUTH, grd%pe_S) call mpp_get_neighbor_pe(grd%domain, EAST, grd%pe_E) call mpp_get_neighbor_pe(grd%domain, WEST, grd%pe_W) !write(stderrunit,'(a,6i4)') 'diamonds, icebergs_init: pe,n,s,e,w =',mpp_pe(),grd%pe_N,grd%pe_S,grd%pe_E,grd%pe_W, NULL_PE !if (verbose) & !write(stderrunit,'(a,i3,a,4i4,a,4f8.2)') 'diamonds, icebergs_init: (',mpp_pe(),') [ij][se]c=', & ! grd%isc,grd%iec,grd%jsc,grd%jec, & ! ' [lon|lat][min|max]=', minval(ice_lon),maxval(ice_lon),minval(ice_lat),maxval(ice_lat) !write(stderrunit,*) 'diamonds, int args = ', mpp_pe(),gni, gnj, layout, axes !write(stderrunit,*) 'diamonds: allocating grid' allocate( grd%lon(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%lon(:,:)=999. allocate( grd%lat(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%lat(:,:)=999. allocate( grd%lonc(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%lon(:,:)=999. allocate( grd%latc(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%lat(:,:)=999. allocate( grd%dx(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%dx(:,:)=0. allocate( grd%dy(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%dy(:,:)=0. allocate( grd%area(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%area(:,:)=0. allocate( grd%msk(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%msk(:,:)=0. allocate( grd%cos(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%cos(:,:)=1. allocate( grd%sin(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%sin(:,:)=0. allocate( grd%calving(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%calving(:,:)=0. allocate( grd%calving_hflx(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%calving_hflx(:,:)=0. allocate( grd%stored_heat(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%stored_heat(:,:)=0. allocate( grd%floating_melt(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%floating_melt(:,:)=0. allocate( grd%berg_melt(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%berg_melt(:,:)=0. allocate( grd%melt_buoy(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%melt_buoy(:,:)=0. allocate( grd%melt_eros(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%melt_eros(:,:)=0. allocate( grd%melt_conv(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%melt_conv(:,:)=0. allocate( grd%bergy_src(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%bergy_src(:,:)=0. allocate( grd%bergy_melt(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%bergy_melt(:,:)=0. allocate( grd%bergy_mass(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%bergy_mass(:,:)=0. allocate( grd%virtual_area(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%virtual_area(:,:)=0. allocate( grd%mass(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%mass(:,:)=0. allocate( grd%mass_on_ocean(grd%isd:grd%ied, grd%jsd:grd%jed, 9) ); grd%mass_on_ocean(:,:,:)=0. allocate( grd%stored_ice(grd%isd:grd%ied, grd%jsd:grd%jed, nclasses) ); grd%stored_ice(:,:,:)=0. allocate( grd%real_calving(grd%isd:grd%ied, grd%jsd:grd%jed, nclasses) ); grd%real_calving(:,:,:)=0. allocate( grd%uo(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%uo(:,:)=0. allocate( grd%vo(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%vo(:,:)=0. allocate( grd%ui(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%ui(:,:)=0. allocate( grd%vi(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%vi(:,:)=0. allocate( grd%ua(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%ua(:,:)=0. allocate( grd%va(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%va(:,:)=0. allocate( grd%ssh(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%ssh(:,:)=0. allocate( grd%sst(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%sst(:,:)=0. allocate( grd%cn(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%cn(:,:)=0. allocate( grd%hi(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%hi(:,:)=0. allocate( grd%tmp(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%tmp(:,:)=0. allocate( grd%tmpc(grd%isc:grd%iec, grd%jsc:grd%jec) ); grd%tmpc(:,:)=0. allocate( bergs%nbergs_calved_by_class(nclasses) ); bergs%nbergs_calved_by_class(:)=0 !write(stderrunit,*) 'diamonds: copying grid' ! Copy data declared on ice model computational domain is=grd%isc; ie=grd%iec; js=grd%jsc; je=grd%jec grd%lon(is:ie,js:je)=ice_lon(:,:) grd%lat(is:ie,js:je)=ice_lat(:,:) grd%area(is:ie,js:je)=ice_area(:,:)*(4.*pi*radius*radius) ! Copy data declared on ice model data domain is=grd%isc-1; ie=grd%iec+1; js=grd%jsc-1; je=grd%jec+1 grd%dx(is:ie,js:je)=ice_dx(:,:) grd%dy(is:ie,js:je)=ice_dy(:,:) grd%msk(is:ie,js:je)=ice_wet(:,:) grd%cos(is:ie,js:je)=cos_rot(:,:) grd%sin(is:ie,js:je)=sin_rot(:,:) call mpp_update_domains(grd%lon, grd%domain, position=CORNER) call mpp_update_domains(grd%lat, grd%domain, position=CORNER) call mpp_update_domains(grd%dy, grd%dx, grd%domain, gridtype=CGRID_NE, flags=SCALAR_PAIR) call mpp_update_domains(grd%area, grd%domain) call mpp_update_domains(grd%msk, grd%domain) call mpp_update_domains(grd%cos, grd%domain, position=CORNER) call mpp_update_domains(grd%sin, grd%domain, position=CORNER) ! Sanitize lon and lat at the SW edges do j=grd%jsc-1,grd%jsd,-1; do i=grd%isd,grd%ied if (grd%lon(i,j).gt.900.) grd%lon(i,j)=grd%lon(i,j+1) if (grd%lat(i,j).gt.900.) grd%lat(i,j)=2.*grd%lat(i,j+1)-grd%lat(i,j+2) enddo; enddo do j=grd%jsd,grd%jed; do i=grd%isd,grd%ied if (grd%lon(i,j).gt.900.) write(stderrunit,*) 'bad lon: ',mpp_pe(),i-grd%isc+1,j-grd%jsc+1 if (grd%lat(i,j).gt.900.) write(stderrunit,*) 'bad lat: ',mpp_pe(),i-grd%isc+1,j-grd%jsc+1 enddo; enddo ! Sanitize lon for the tile (need continuous longitudes within one tile) j=grd%jsc; do i=grd%isc+1,grd%ied minl=grd%lon(i-1,j)-180. grd%lon(i,j)=modulo(grd%lon(i,j)-minl,360.)+minl enddo j=grd%jsc; do i=grd%isc-1,grd%isd,-1 minl=grd%lon(i+1,j)-180. grd%lon(i,j)=modulo(grd%lon(i,j)-minl,360.)+minl enddo do j=grd%jsc+1,grd%jed; do i=grd%isd,grd%ied minl=grd%lon(i,j-1)-180. grd%lon(i,j)=modulo(grd%lon(i,j)-minl,360.)+minl enddo; enddo do j=grd%jsc-1,grd%jsd,-1; do i=grd%isd,grd%ied minl=grd%lon(i,j+1)-180. grd%lon(i,j)=modulo(grd%lon(i,j)-minl,360.)+minl enddo; enddo ! lonc, latc used for searches do j=grd%jsd+1,grd%jed; do i=grd%isd+1,grd%ied grd%lonc(i,j)=0.25*( (grd%lon(i,j)+grd%lon(i-1,j-1)) & +(grd%lon(i-1,j)+grd%lon(i,j-1)) ) grd%latc(i,j)=0.25*( (grd%lat(i,j)+grd%lat(i-1,j-1)) & +(grd%lat(i-1,j)+grd%lat(i,j-1)) ) enddo; enddo if (debug) then write(stderrunit,'(a,i3,a,4i4,a,4f8.2)') 'diamonds, icebergs_init: (',mpp_pe(),') [ij][se]c=', & grd%isc,grd%iec,grd%jsc,grd%jec, & ' [lon|lat][min|max]=', minval(grd%lon),maxval(grd%lon),minval(grd%lat),maxval(grd%lat) endif !if (mpp_pe().eq.3) then ! write(stderrunit,'(a3,32i7)') 'Lon',(i,i=grd%isd,grd%ied) ! do j=grd%jed,grd%jsd,-1 ! write(stderrunit,'(i3,32f7.1)') j,(grd%lon(i,j),i=grd%isd,grd%ied) ! enddo ! write(stderrunit,'(a3,32i7)') 'Lat',(i,i=grd%isd,grd%ied) ! do j=grd%jed,grd%jsd,-1 ! write(stderrunit,'(i3,32f7.1)') j,(grd%lat(i,j),i=grd%isd,grd%ied) ! enddo ! write(stderrunit,'(a3,32i7)') 'Msk',(i,i=grd%isd,grd%ied) ! do j=grd%jed,grd%jsd,-1 ! write(stderrunit,'(i3,32f7.1)') j,(grd%msk(i,j),i=grd%isd,grd%ied) ! enddo !endif ! Parameters bergs%dt=dt bergs%traj_sample_hrs=traj_sample_hrs bergs%verbose_hrs=verbose_hrs bergs%grd%halo=halo bergs%rho_bergs=rho_bergs bergs%LoW_ratio=LoW_ratio bergs%use_operator_splitting=use_operator_splitting bergs%bergy_bit_erosion_fraction=bergy_bit_erosion_fraction bergs%sicn_shift=sicn_shift bergs%passive_mode=passive_mode bergs%time_average_weight=time_average_weight bergs%speed_limit=speed_limit bergs%add_weight_to_ocean=add_weight_to_ocean allocate( bergs%initial_mass(nclasses) ); bergs%initial_mass(:)=initial_mass(:) allocate( bergs%distribution(nclasses) ); bergs%distribution(:)=distribution(:) allocate( bergs%mass_scaling(nclasses) ); bergs%mass_scaling(:)=mass_scaling(:) allocate( bergs%initial_thickness(nclasses) ); bergs%initial_thickness(:)=initial_thickness(:) allocate( bergs%initial_width(nclasses) ) allocate( bergs%initial_length(nclasses) ) bergs%initial_width(:)=sqrt(initial_mass(:)/(LoW_ratio*rho_bergs*initial_thickness(:))) bergs%initial_length(:)=LoW_ratio*bergs%initial_width(:) call mpp_clock_end(bergs%clock_ini) call mpp_clock_begin(bergs%clock_ior) call read_restart_bergs(bergs) call bergs_chksum(bergs, 'read_restart bergs') call read_restart_calving(bergs) call mpp_clock_end(bergs%clock_ior) call mpp_clock_begin(bergs%clock_ini) if (really_debug) call print_bergs(stderrunit,bergs,'icebergs_init, initial status') ! Diagnostics id_class = diag_axis_init('mass_class', initial_mass, 'kg','Z', 'iceberg mass') axes3d(1:2)=axes axes3d(3)=id_class grd%id_calving=register_diag_field('icebergs', 'calving', axes, Time, & 'Incoming Calving mass rate', 'kg/s') grd%id_calving_hflx_in=register_diag_field('icebergs', 'calving_hflx_in', axes, Time, & 'Incoming Calving heat flux', 'J/s') grd%id_accum=register_diag_field('icebergs', 'accum_calving', axes, Time, & 'Accumulated calving mass rate', 'kg/s') grd%id_unused=register_diag_field('icebergs', 'unused_calving', axes, Time, & 'Unused calving mass rate', 'kg/s') grd%id_floating_melt=register_diag_field('icebergs', 'melt', axes, Time, & 'Melt rate of icebergs + bits', 'kg/(m^2*s)') grd%id_berg_melt=register_diag_field('icebergs', 'berg_melt', axes, Time, & 'Melt rate of icebergs', 'kg/(m^2*s)') grd%id_melt_buoy=register_diag_field('icebergs', 'melt_buoy', axes, Time, & 'Buoyancy component of iceberg melt rate', 'kg/(m^2*s)') grd%id_melt_eros=register_diag_field('icebergs', 'melt_eros', axes, Time, & 'Erosion component of iceberg melt rate', 'kg/(m^2*s)') grd%id_melt_conv=register_diag_field('icebergs', 'melt_conv', axes, Time, & 'Convective component of iceberg melt rate', 'kg/(m^2*s)') grd%id_bergy_src=register_diag_field('icebergs', 'bergy_src', axes, Time, & 'Mass source of bergy bits', 'kg/(m^2*s)') grd%id_bergy_melt=register_diag_field('icebergs', 'bergy_melt', axes, Time, & 'Melt rate of bergy bits', 'kg/(m^2*s)') grd%id_bergy_mass=register_diag_field('icebergs', 'bergy_mass', axes, Time, & 'Bergy bit density field', 'kg/(m^2)') grd%id_virtual_area=register_diag_field('icebergs', 'virtual_area', axes, Time, & 'Virtual coverage by icebergs', 'm^2') grd%id_mass=register_diag_field('icebergs', 'mass', axes, Time, & 'Iceberg density field', 'kg/(m^2)') grd%id_stored_ice=register_diag_field('icebergs', 'stored_ice', axes3d, Time, & 'Accumulated ice mass by class', 'kg') grd%id_real_calving=register_diag_field('icebergs', 'real_calving', axes3d, Time, & 'Calving into iceberg class', 'kg/s') grd%id_uo=register_diag_field('icebergs', 'uo', axes, Time, & 'Ocean zonal component of velocity', 'm s^-1') grd%id_vo=register_diag_field('icebergs', 'vo', axes, Time, & 'Ocean meridional component of velocity', 'm s^-1') grd%id_ui=register_diag_field('icebergs', 'ui', axes, Time, & 'Ice zonal component of velocity', 'm s^-1') grd%id_vi=register_diag_field('icebergs', 'vi', axes, Time, & 'Ice meridional component of velocity', 'm s^-1') grd%id_ua=register_diag_field('icebergs', 'ua', axes, Time, & 'Atmos zonal component of velocity', 'm s^-1') grd%id_va=register_diag_field('icebergs', 'va', axes, Time, & 'Atmos meridional component of velocity', 'm s^-1') grd%id_sst=register_diag_field('icebergs', 'sst', axes, Time, & 'Sea surface temperature', 'degrees_C') grd%id_cn=register_diag_field('icebergs', 'cn', axes, Time, & 'Sea ice concentration', '(fraction)') grd%id_hi=register_diag_field('icebergs', 'hi', axes, Time, & 'Sea ice thickness', 'm') ! Static fields id_class=register_static_field('icebergs', 'lon', axes, & 'longitude (corners)', 'degrees_E') if (id_class>0) lerr=send_data(id_class, grd%lon(grd%isc:grd%iec,grd%jsc:grd%jec)) id_class=register_static_field('icebergs', 'lat', axes, & 'latitude (corners)', 'degrees_N') if (id_class>0) lerr=send_data(id_class, grd%lat(grd%isc:grd%iec,grd%jsc:grd%jec)) id_class=register_static_field('icebergs', 'area', axes, & 'cell area', 'm^2') if (id_class>0) lerr=send_data(id_class, grd%area(grd%isc:grd%iec,grd%jsc:grd%jec)) id_class=register_static_field('icebergs', 'mask', axes, & 'wet point mask', 'none') if (id_class>0) lerr=send_data(id_class, grd%msk(grd%isc:grd%iec,grd%jsc:grd%jec)) if (debug) then call grd_chksum2(grd, grd%lon, 'init lon') call grd_chksum2(grd, grd%lat, 'init lat') call grd_chksum2(grd, grd%area, 'init area') call grd_chksum2(grd, grd%msk, 'init msk') call grd_chksum2(grd, grd%cos, 'init cos') call grd_chksum2(grd, grd%sin, 'init sin') endif !write(stderrunit,*) 'diamonds: done' call mpp_clock_end(bergs%clock_ini) call mpp_clock_end(bergs%clock) end subroutine icebergs_init ! ############################################################################## subroutine read_restart_bergs(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables integer :: k, ierr, ncid, dimid, nbergs_in_file integer :: lonid, latid, uvelid, vvelid, ineid, jneid integer :: massid, thicknessid, widthid, lengthid integer :: start_lonid, start_latid, start_yearid, start_dayid, start_massid integer :: scaling_id, mass_of_bits_id, heat_density_id logical :: lres, found_restart, multiPErestart=.FALSE. real :: lon0, lon1, lat0, lat1 character(len=30) :: filename type(icebergs_gridded), pointer :: grd type(iceberg) :: localberg ! NOT a pointer but an actual local variable integer :: stderrunit stderrunit=stderr() ! For convenience grd=>bergs%grd ! Find a restart file multiPErestart=.FALSE. do filename='INPUT/icebergs.res.nc'; inquire(file=filename,exist=found_restart) if (found_restart) exit write(filename(1:27),'("INPUT/icebergs.res.nc.",I4.4)') mpp_pe() inquire(file=filename,exist=found_restart) if (found_restart) multiPErestart=.TRUE. if (found_restart) exit write(filename(1:21),'("icebergs.res.nc.",I4.4)') mpp_pe() inquire(file=filename,exist=found_restart) if (found_restart) multiPErestart=.TRUE. if (found_restart) exit filename='icebergs.res.nc'; inquire(file=filename,exist=found_restart) if (found_restart) exit if (verbose.and.mpp_pe()==mpp_root_pe()) write(*,'(a)') 'diamonds, read_restart_bergs: no restart file found' ! return ! leave s/r if no restart found multiPErestart=.TRUE. ! This is to force sanity checking in a mulit-PE mode if no file was found on this PE exit enddo if (found_restart) then ! only do the following if a file was found if (verbose.and.mpp_pe()==mpp_root_pe()) write(*,'(2a)') 'diamonds, read_restart_bergs: found restart file = ',filename ierr=nf_open(filename, NF_NOWRITE, ncid) if (ierr .ne. NF_NOERR) write(stderrunit,*) 'diamonds, read_restart_bergs: nf_open failed' ierr=nf_inq_unlimdim(ncid, dimid) if (ierr .ne. NF_NOERR) write(stderrunit,*) 'diamonds, read_restart_bergs: nf_inq_unlimdim failed' ierr=nf_inq_dimlen(ncid, dimid, nbergs_in_file) if (ierr .ne. NF_NOERR) write(stderrunit,*) 'diamonds, read_restart_bergs: nf_inq_dimlen failed' !write(stderrunit,*) 'diamonds, read_restart_bergs: nbergs in file =', nbergs_in_file lonid=inq_var(ncid, 'lon') latid=inq_var(ncid, 'lat') uvelid=inq_var(ncid, 'uvel') vvelid=inq_var(ncid, 'vvel') massid=inq_var(ncid, 'mass') thicknessid=inq_var(ncid, 'thickness') widthid=inq_var(ncid, 'width') lengthid=inq_var(ncid, 'length') start_lonid=inq_var(ncid, 'start_lon') start_latid=inq_var(ncid, 'start_lat') start_yearid=inq_var(ncid, 'start_year') start_dayid=inq_var(ncid, 'start_day') start_massid=inq_var(ncid, 'start_mass') scaling_id=inq_var(ncid, 'mass_scaling') mass_of_bits_id=inq_var(ncid, 'mass_of_bits',unsafe=.true.) heat_density_id=inq_var(ncid, 'heat_density',unsafe=.true.) ineid=inq_var(ncid, 'ine',unsafe=.true.) jneid=inq_var(ncid, 'jne',unsafe=.true.) ! Find approx outer bounds for tile lon0=minval( grd%lon(grd%isc-1:grd%iec,grd%jsc-1:grd%jec) ) lon1=maxval( grd%lon(grd%isc-1:grd%iec,grd%jsc-1:grd%jec) ) lat0=minval( grd%lat(grd%isc-1:grd%iec,grd%jsc-1:grd%jec) ) lat1=maxval( grd%lat(grd%isc-1:grd%iec,grd%jsc-1:grd%jec) ) !write(stderrunit,'(a,i3,a,4f10.3)') 'diamonds, read_restart_bergs: (',mpp_pe(),') ',lon0,lon1,lat0,lat1 do k=1, nbergs_in_file !write(stderrunit,*) 'diamonds, read_restart_bergs: reading berg ',k localberg%lon=get_double(ncid, lonid, k) localberg%lat=get_double(ncid, latid, k) if (ineid>0 .and. jneid>0 .and. .not. ignore_ij_restart) then ! read i,j position and avoid the "find" step localberg%ine=get_int(ncid, ineid, k) localberg%jne=get_int(ncid, jneid, k) if ( localberg%ine>=grd%isc .and. localberg%ine<=grd%iec .and. & localberg%jne>=grd%jsc .and.localberg%jne<=grd%jec ) then lres=.true. else lres=.false. endif else ! i,j are not available from the file so we search the grid to find out if we reside on this PE if (use_slow_find) then lres=find_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne) else lres=find_cell_by_search(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne) endif endif if (really_debug) then write(stderrunit,'(a,i,a,2f,a,i)') 'diamonds, read_restart_bergs: berg ',k,' is at ',localberg%lon,localberg%lat,' on PE ',mpp_pe() write(stderrunit,*) 'diamonds, read_restart_bergs: lres = ',lres endif if (lres) then ! true if we reside on this PE grid localberg%uvel=get_double(ncid, uvelid, k) localberg%vvel=get_double(ncid, vvelid, k) localberg%mass=get_double(ncid, massid, k) localberg%thickness=get_double(ncid, thicknessid, k) localberg%width=get_double(ncid, widthid, k) localberg%length=get_double(ncid, lengthid, k) localberg%start_lon=get_double(ncid, start_lonid, k) localberg%start_lat=get_double(ncid, start_latid, k) localberg%start_year=get_int(ncid, start_yearid, k) localberg%start_day=get_double(ncid, start_dayid, k) localberg%start_mass=get_double(ncid, start_massid, k) localberg%mass_scaling=get_double(ncid, scaling_id, k) if (mass_of_bits_id>0) then ! Allow reading of older restart with no bergy bits localberg%mass_of_bits=get_double(ncid, mass_of_bits_id, k) else localberg%mass_of_bits=0. endif if (heat_density_id>0) then ! Allow reading of older restart with no heat content localberg%heat_density=get_double(ncid, heat_density_id, k) else localberg%heat_density=0. endif if (really_debug) lres=is_point_in_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne, explain=.true.) lres=pos_within_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne, localberg%xi, localberg%yj) !call add_new_berg_to_list(bergs%first, localberg, quick=.true.) call add_new_berg_to_list(bergs%first, localberg) if (really_debug) call print_berg(stderrunit, bergs%first, 'read_restart_bergs, add_new_berg_to_list') elseif (multiPErestart) then call error_mesg('diamonds, read_restart_bergs', 'berg in PE file was not on PE!', FATAL) endif enddo else ! if no restart file was read on this PE nbergs_in_file=0 endif ! if (found_restart) ! Sanity check k=count_bergs(bergs) if (verbose) write(*,'(2(a,i))') 'diamonds, read_restart_bergs: # bergs =',k,' on PE',mpp_pe() if (multiPErestart) call mpp_sum(nbergs_in_file) ! In case PE 0 didn't open a file call mpp_sum(k) bergs%nbergs_start=k if (mpp_pe().eq.mpp_root_pe()) then write(*,'(a,i,a,i,a)') 'diamonds, read_restart_bergs: there were',nbergs_in_file,' bergs in the restart file and', & k,' bergs have been read' endif if (k.ne.nbergs_in_file) call error_mesg('diamonds, read_restart_bergs', 'wrong number of bergs read!', FATAL) if (.not. found_restart .and. bergs%nbergs_start==0 .and. generate_test_icebergs) call genereate_bergs(bergs) bergs%floating_mass_start=sum_mass(bergs%first) call mpp_sum( bergs%floating_mass_start ) bergs%icebergs_mass_start=sum_mass(bergs%first,justbergs=.true.) call mpp_sum( bergs%icebergs_mass_start ) bergs%bergy_mass_start=sum_mass(bergs%first,justbits=.true.) call mpp_sum( bergs%bergy_mass_start ) if (mpp_pe().eq.mpp_root_pe().and.verbose) write(*,'(a)') 'diamonds, read_restart_bergs: completed' contains subroutine genereate_bergs(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables integer :: i,j type(iceberg) :: localberg ! NOT a pointer but an actual local variable ! For convenience grd=>bergs%grd do j=grd%jsc,grd%jec; do i=grd%isc,grd%iec if (grd%msk(i,j)>0. .and. abs(grd%latc(i,j))>60.) then localberg%xi=0.5 localberg%yj=0.5 localberg%ine=i localberg%jne=j localberg%lon=bilin(grd, grd%lon, i, j, localberg%xi, localberg%yj) localberg%lat=bilin(grd, grd%lat, i, j, localberg%xi, localberg%yj) localberg%mass=bergs%initial_mass(1) localberg%thickness=bergs%initial_thickness(1) localberg%width=bergs%initial_width(1) localberg%length=bergs%initial_length(1) localberg%start_lon=localberg%lon localberg%start_lat=localberg%lat localberg%start_year=0 localberg%start_day=1 localberg%start_mass=localberg%mass localberg%mass_scaling=bergs%mass_scaling(1) localberg%mass_of_bits=0. localberg%heat_density=0. localberg%uvel=1. localberg%vvel=0. call add_new_berg_to_list(bergs%first, localberg) localberg%uvel=-1. localberg%vvel=0. call add_new_berg_to_list(bergs%first, localberg) localberg%uvel=0. localberg%vvel=1. call add_new_berg_to_list(bergs%first, localberg) localberg%uvel=0. localberg%vvel=-1. call add_new_berg_to_list(bergs%first, localberg) endif enddo; enddo bergs%nbergs_start=count_bergs(bergs) call mpp_sum(bergs%nbergs_start) if (mpp_pe().eq.mpp_root_pe()) & write(*,'(a,i,a)') 'diamonds, generate_bergs: ',bergs%nbergs_start,' were generated' end subroutine genereate_bergs end subroutine read_restart_bergs ! ############################################################################## subroutine read_restart_calving(bergs) use random_numbers_mod, only: initializeRandomNumberStream, getRandomNumbers, randomNumberStream ! Arguments type(icebergs), pointer :: bergs ! Local variables integer :: k,i character(len=30) :: filename type(icebergs_gridded), pointer :: grd real, allocatable, dimension(:,:) :: randnum type(randomNumberStream) :: rns ! For convenience grd=>bergs%grd ! Read stored ice filename='INPUT/calving.res.nc' if (file_exist(filename)) then if (verbose.and.mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') & 'diamonds, read_restart_calving: reading ',filename call read_data(filename, 'stored_ice', grd%stored_ice, grd%domain) if (field_exist(filename, 'stored_heat')) then if (verbose.and.mpp_pe().eq.mpp_root_pe()) write(*,'(a)') & 'diamonds, read_restart_calving: reading stored_heat from restart file.' call read_data(filename, 'stored_heat', grd%stored_heat, grd%domain) else if (verbose.and.mpp_pe().eq.mpp_root_pe()) write(*,'(a)') & 'diamonds, read_restart_calving: stored_heat WAS NOT FOUND in the file. Setting to 0.' grd%stored_heat(:,:)=0. endif bergs%restarted=.true. else if (verbose.and.mpp_pe().eq.mpp_root_pe()) write(*,'(a)') & 'diamonds, read_restart_calving: initializing stored ice to random numbers' allocate(randnum(grd%jsc:grd%jec, nclasses)) do i=grd%isc, grd%iec rns=initializeRandomNumberStream(i) call getRandomNumbers(rns,randnum) do k=1, nclasses grd%stored_ice(i,grd%jsc:grd%jec,k)=randnum(:,k) & *grd%msk(i,grd%jsc:grd%jec) & *bergs%initial_mass(k)*bergs%mass_scaling(k) enddo enddo deallocate(randnum) endif call grd_chksum3(bergs%grd, bergs%grd%stored_ice, 'read_restart_calving, stored_ice') call grd_chksum2(bergs%grd, bergs%grd%stored_heat, 'read_restart_calving, stored_heat') bergs%stored_start=sum( grd%stored_ice(grd%isc:grd%iec,grd%jsc:grd%jec,:) ) call mpp_sum( bergs%stored_start ) bergs%stored_heat_start=sum( grd%stored_heat(grd%isc:grd%iec,grd%jsc:grd%jec) ) call mpp_sum( bergs%stored_heat_start ) bergs%floating_heat_start=sum_heat(bergs%first) call mpp_sum( bergs%floating_heat_start ) end subroutine read_restart_calving ! ############################################################################## subroutine add_new_berg_to_list(first, bergvals, quick) ! Arguments type(iceberg), pointer :: first type(iceberg), intent(in) :: bergvals logical, intent(in), optional :: quick ! Local variables type(iceberg), pointer :: new=>NULL() new=>NULL() call create_iceberg(new, bergvals) if (present(quick).and. quick) then call insert_berg_into_list(first, new, quick=.true.) else call insert_berg_into_list(first, new) endif !Clear new new=>NULL() end subroutine add_new_berg_to_list ! ############################################################################## subroutine count_out_of_order(bergs,label) ! Arguments type(icebergs), pointer :: bergs character(len=*) :: label ! Local variables type(iceberg), pointer :: this, next integer :: i, icnt1, icnt2 icnt1=0 this=>bergs%first next=>NULL() if (associated(this)) then if (associated(this%next)) next=>this%next endif do while (associated(next)) if (.not. inorder(this,next)) icnt1=icnt1+1 this=>next next=>next%next enddo call mpp_sum(icnt1) i=0; icnt2=0 this=>bergs%first do while (associated(this)) i=1 if (this%inebergs%grd%iec .or. & this%jnebergs%grd%jec) icnt2=icnt2+1 this=>this%next if (i>1.and..not.associated(this%prev)) then call error_mesg('diamonds, count_out_of_order', 'Pointer %prev is unassociated. This should not happen!', FATAL) endif enddo call mpp_sum(icnt2) if ((debug.or.icnt1.ne.0).and.mpp_pe().eq.mpp_root_pe()) then write(*,'(a,2(x,a,i6),x,a)') 'diamonds, count_out_of_order:', & '# out of order=', icnt1,'# in halo=',icnt2,label endif call check_for_duplicates(bergs,label) end subroutine count_out_of_order ! ############################################################################## subroutine check_for_duplicates(bergs,label) ! Arguments type(icebergs), pointer :: bergs character(len=*) :: label ! Local variables type(iceberg), pointer :: this1, next1, this2, next2 integer :: icnt_id, icnt_same icnt_id=0 icnt_same=0 this1=>bergs%first do while (associated(this1)) this2=>this1%next do while (associated(this2)) if (sameid(this1,this2)) icnt_id=icnt_id+1 if (sameberg(this1,this2)) icnt_same=icnt_same+1 this2=>this2%next enddo this1=>this1%next enddo call mpp_sum(icnt_id) call mpp_sum(icnt_same) if ((debug.or.icnt_id>0.or.icnt_same>0).and.mpp_pe().eq.mpp_root_pe()) then write(*,'(a,2(x,a,i9),x,a)') 'diamonds, check_for_duplicates:', & '# with same id=', icnt_id,'# identical bergs=',icnt_same,label endif end subroutine check_for_duplicates ! ############################################################################## subroutine insert_berg_into_list(first, newberg, quick) ! Arguments type(iceberg), pointer :: first, newberg logical, intent(in), optional :: quick ! Local variables type(iceberg), pointer :: this, prev if (associated(first)) then if (.not. parallel_reprod .or. (present(quick).and.quick)) then newberg%next=>first first%prev=>newberg first=>newberg else if (inorder(newberg,first)) then ! Insert at front of list newberg%next=>first first%prev=>newberg first=>newberg else this=>first prev=>NULL() do while( associated(this) ) if (inorder(newberg,this) ) then exit endif prev=>this this=>this%next enddo prev%next=>newberg newberg%prev=>prev if (associated(this)) this%prev=>newberg newberg%next=>this endif endif else ! list is empty so create it first=>newberg endif end subroutine insert_berg_into_list ! ############################################################################## logical function inorder(berg1, berg2) ! Arguments type(iceberg), pointer :: berg1, berg2 ! Local variables if (berg1%start_yearberg2%start_year) then inorder=.false. return endif if (berg1%start_dayberg2%start_day) then inorder=.false. return endif if (berg1%start_massberg2%start_mass) then inorder=.false. return endif if (berg1%start_lonberg2%start_lon) then inorder=.false. return endif if (berg1%start_latberg2%start_lat) then inorder=.false. return endif inorder=.true. ! passing the above tests mean the bergs 1 and 2 are identical? end function inorder ! ############################################################################## real function time_hash(berg) ! Arguments type(iceberg), pointer :: berg time_hash=berg%start_day+366.*float(berg%start_year) end function time_hash ! ############################################################################## real function pos_hash(berg) ! Arguments type(iceberg), pointer :: berg pos_hash=berg%start_lon+360.*(berg%start_lat+90.) end function pos_hash ! ############################################################################## logical function sameid(berg1, berg2) ! Arguments type(iceberg), pointer :: berg1, berg2 ! Local variables sameid=.false. if (berg1%start_year.ne.berg2%start_year) return if (berg1%start_day.ne.berg2%start_day) return if (berg1%start_mass.ne.berg2%start_mass) return if (berg1%start_lon.ne.berg2%start_lon) return if (berg1%start_lat.ne.berg2%start_lat) return sameid=.true. ! passing the above tests means that bergs 1 and 2 have the same id end function sameid ! ############################################################################## logical function sameberg(berg1, berg2) ! Arguments type(iceberg), pointer :: berg1, berg2 ! Local variables sameberg=.false. if (.not. sameid(berg1, berg2)) return if (berg1%lon.ne.berg2%lon) return if (berg1%lat.ne.berg2%lat) return if (berg1%mass.ne.berg2%mass) return if (berg1%uvel.ne.berg2%uvel) return if (berg1%vvel.ne.berg2%vvel) return if (berg1%thickness.ne.berg2%thickness) return if (berg1%width.ne.berg2%width) return if (berg1%length.ne.berg2%length) return sameberg=.true. ! passing the above tests mean that bergs 1 and 2 are identical end function sameberg ! ############################################################################## subroutine create_iceberg(berg, bergvals) ! Arguments type(iceberg), pointer :: berg type(iceberg), intent(in) :: bergvals ! Local variables integer :: stderrunit stderrunit=stderr() if (associated(berg)) then write(stderrunit,*) 'diamonds, create_iceberg: berg already associated!!!!',mpp_pe() call error_mesg('diamonds, create_iceberg', 'berg already associated. This should not happen!', FATAL) endif allocate(berg) berg=bergvals berg%prev=>NULL() berg%next=>NULL() end subroutine create_iceberg ! ############################################################################## subroutine delete_iceberg_from_list(first, berg) ! Arguments type(iceberg), pointer :: first, berg ! Local variables ! Connect neighbors to each other if (associated(berg%prev)) then berg%prev%next=>berg%next else first=>berg%next endif if (associated(berg%next)) berg%next%prev=>berg%prev ! Bye-bye berg call destroy_iceberg(berg) end subroutine delete_iceberg_from_list ! ############################################################################## subroutine destroy_iceberg(berg) ! Arguments type(iceberg), pointer :: berg ! Local variables ! Bye-bye berg deallocate(berg) end subroutine destroy_iceberg ! ############################################################################## subroutine print_berg(iochan, berg, label) ! Arguments integer, intent(in) :: iochan type(iceberg), pointer :: berg character(len=*) :: label ! Local variables write(iochan,'("diamonds, print_berg: ",a," pe=(",i3,") start lon,lat,yr,day,mass=",2f10.4,i5,f7.2,es12.4)') & label, mpp_pe(), berg%start_lon, berg%start_lat, & berg%start_year, berg%start_day, berg%start_mass write(iochan,'("diamonds, print_berg: ",a," pe=(",i3,a,2i5,3(a,2f10.4),a,2l2)') & label, mpp_pe(), ') i,j=',berg%ine, berg%jne, & ' xi,yj=', berg%xi, berg%yj, & ' lon,lat=', berg%lon, berg%lat, & ' u,v=', berg%uvel, berg%vvel, & ' p,n=', associated(berg%prev), associated(berg%next) write(iochan,'("diamonds, print_berg: ",a," pe=(",i3,") ",6(a,2f10.4))') & label, mpp_pe(), 'uo,vo=', berg%uo, berg%vo, 'ua,va=', berg%ua, berg%va, 'ui,vi=', berg%ui, berg%vi end subroutine print_berg ! ############################################################################## subroutine print_bergs(iochan, bergs, label) ! Arguments integer, intent(in) :: iochan type(icebergs), pointer :: bergs character(len=*) :: label ! Local variables integer :: nbergs, nnbergs type(iceberg), pointer :: this this=>bergs%first do while(associated(this)) call print_berg(iochan, this, label) this=>this%next enddo nbergs=count_bergs(bergs) nnbergs=nbergs call mpp_sum(nnbergs) if (nbergs.gt.0) write(iochan,'("diamonds, ",a," there are",i5," bergs out of",i6," on PE ",i4)') label, nbergs, nnbergs, mpp_pe() end subroutine print_bergs ! ############################################################################## integer function count_bergs(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables type(iceberg), pointer :: this count_bergs=0 this=>bergs%first do while(associated(this)) count_bergs=count_bergs+1 this=>this%next enddo end function count_bergs ! ############################################################################## subroutine record_posn(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables type(xyt) :: posn type(iceberg), pointer :: this this=>bergs%first do while (associated(this)) posn%lon=this%lon posn%lat=this%lat posn%year=bergs%current_year posn%day=bergs%current_yearday posn%uvel=this%uvel posn%vvel=this%vvel posn%mass=this%mass posn%mass_of_bits=this%mass_of_bits posn%heat_density=this%heat_density posn%thickness=this%thickness posn%width=this%width posn%length=this%length posn%uo=this%uo posn%vo=this%vo posn%ui=this%ui posn%vi=this%vi posn%ua=this%ua posn%va=this%va posn%ssh_x=this%ssh_x posn%ssh_y=this%ssh_y posn%sst=this%sst posn%cn=this%cn posn%hi=this%hi call push_posn(this%trajectory, posn) this=>this%next enddo end subroutine record_posn ! ############################################################################## subroutine push_posn(trajectory, posn_vals) ! Arguments type(xyt), pointer :: trajectory type(xyt) :: posn_vals ! Local variables type(xyt), pointer :: new_posn allocate(new_posn) new_posn=posn_vals new_posn%next=>trajectory trajectory=>new_posn end subroutine push_posn ! ############################################################################## subroutine move_trajectory(bergs, berg) ! Arguments type(icebergs), pointer :: bergs type(iceberg), pointer :: berg ! Local variables type(xyt), pointer :: next, last type(xyt) :: vals ! If the trajectory is empty, ignore it if (.not.associated(berg%trajectory)) return ! Push identifying info into first posn (note reverse order due to stack) vals%lon=berg%start_lon vals%lat=berg%start_lat vals%year=berg%start_year vals%day=berg%start_day vals%mass=berg%start_mass call push_posn(berg%trajectory, vals) vals%lon=0. vals%lat=99. vals%year=0 vals%day=0. vals%mass=berg%mass_scaling call push_posn(berg%trajectory, vals) ! Find end of berg trajectory and point it to start of existing trajectories next=>berg%trajectory do while (associated(next)) last=>next next=>next%next enddo last%next=>bergs%trajectories bergs%trajectories=>berg%trajectory berg%trajectory=>NULL() end subroutine move_trajectory ! ############################################################################## logical function find_cell_by_search(grd, x, y, i, j) ! Arguments type(icebergs_gridded), pointer :: grd real, intent(in) :: x, y integer, intent(inout) :: i, j ! Local variables integer :: is,ie,js,je,di,dj,io,jo,icnt real :: d0,d1,d2,d3,d4,d5,d6,d7,d8,dmin logical :: explain=.false. 911 continue find_cell_by_search=.false. is=grd%isc; ie=grd%iec; js=grd%jsc; je=grd%jec ! Start at nearest corner d1=dcost(x,y,grd%lonc(is+1,js+1),grd%latc(is+1,js+1)) d2=dcost(x,y,grd%lonc(ie-1,js+1),grd%latc(ie-1,js+1)) d3=dcost(x,y,grd%lonc(ie-1,je-1),grd%latc(ie-1,je-1)) d4=dcost(x,y,grd%lonc(is+1,je-1),grd%latc(is+1,je-1)) dmin=min(d1,d2,d3,d4) if (d1==dmin) then; i=is+1; j=js+1 elseif (d2==dmin) then; i=ie-1; j=js+1 elseif (d3==dmin) then; i=ie-1; j=je-1 elseif (d4==dmin) then; i=is+1; j=je-1 else call error_mesg('diamonds, find_cell_by_search:', 'This should never EVER happen! (1)', FATAL) endif if (explain) then write(0,'(i3,a,2i4,f)') mpp_pe(),'Initial corner i-is,j-js,cost=',i-is,j-js,dmin write(0,'(i3,a,3f)') mpp_pe(),'cost ',d4,d3 write(0,'(i3,a,3f)') mpp_pe(),'cost ',d1,d2 endif if (is_point_in_cell(grd, x, y, i, j)) then find_cell_by_search=.true. return endif do icnt=1, 1*(ie-is+je-js) io=i; jo=j d0=dcost(x,y,grd%lonc(io,jo),grd%latc(io,jo)) d1=dcost(x,y,grd%lonc(io,jo+1),grd%latc(io,jo+1)) d2=dcost(x,y,grd%lonc(io-1,jo+1),grd%latc(io-1,jo+1)) d3=dcost(x,y,grd%lonc(io-1,jo),grd%latc(io-1,jo)) d4=dcost(x,y,grd%lonc(io-1,jo-1),grd%latc(io-1,jo-1)) d5=dcost(x,y,grd%lonc(io,jo-1),grd%latc(io,jo-1)) d6=dcost(x,y,grd%lonc(io+1,jo-1),grd%latc(io+1,jo-1)) d7=dcost(x,y,grd%lonc(io+1,jo),grd%latc(io+1,jo)) d8=dcost(x,y,grd%lonc(io+1,jo+1),grd%latc(io+1,jo+1)) ! dmin=min(d0,d1,d3,d5,d7) dmin=min(d0,d1,d2,d3,d4,d5,d6,d7,d8) if (d0==dmin) then; di=0; dj=0 elseif (d2==dmin) then; di=-1; dj=1 elseif (d4==dmin) then; di=-1; dj=-1 elseif (d6==dmin) then; di=1; dj=-1 elseif (d8==dmin) then; di=1; dj=1 elseif (d1==dmin) then; di=0; dj=1 elseif (d3==dmin) then; di=-1; dj=0 elseif (d5==dmin) then; di=0; dj=-1 elseif (d7==dmin) then; di=1; dj=0 else call error_mesg('diamonds, find_cell_by_search:', 'This should never EVER happen!', FATAL) endif i=min(ie, max(is, io+di)) j=min(je, max(js, jo+dj)) if (explain) then write(0,'(i3,a,2i4,f,a,2i4,a,2i4)') mpp_pe(),'Current position i,j,cost=',i,j,dmin,' di,dj=',di,dj,' old io,jo=',io,jo write(0,'(i3,a,3f)') mpp_pe(),'cost ',d2,d1,d8 write(0,'(i3,a,3f)') mpp_pe(),'cost ',d3,d0,d7 write(0,'(i3,a,3f)') mpp_pe(),'cost ',d4,d5,d6 endif if (is_point_in_cell(grd, x, y, i, j)) then find_cell_by_search=.true. return endif if ((i==io.and.j==jo) & .and. .not.find_better_min(grd, x, y, 3, i, j) & ) then ! Stagnated find_cell_by_search=find_cell_loc(grd, x, y, is, ie, js, je, 1, i, j) if (.not. find_cell_by_search) find_cell_by_search=find_cell_loc(grd, x, y, is, ie, js, je, 3, i, j) i=min(ie, max(is, i)) j=min(je, max(js, j)) if (is_point_in_cell(grd, x, y, i, j)) then find_cell_by_search=.true. else ! find_cell_by_search=find_cell(grd, x, y, io, jo) ! if (find_cell_by_search) then ! if (explain) then ! call print_fld(grd, grd%lat, 'Lat') ! call print_fld(grd, grd%lon, 'Lat') ! do j=grd%jsd, grd%jed; do i=grd%isd, grd%ied ! grd%tmp(i,j) = dcost(x,y,grd%lonc(i,j),grd%latc(i,j)) ! enddo; enddo ! call print_fld(grd, grd%tmp, 'Cost') ! stop 'Avoid recursing' ! endif ! write(0,'(i3,a,2i5,a,2i3,a,2f8.3)') mpp_pe(),'diamonds, find_cell_by_search: false negative io,jo=',io,jo,' di,dj=',di,dj,' targ=',x,y ! explain=.true.; goto 911 ! endif endif return endif enddo find_cell_by_search=find_cell(grd, x, y, i, j) if (find_cell_by_search) then write(0,'(i3,a,3f)') mpp_pe(),'cost ',d2,d1,d8 write(0,'(i3,a,3f)') mpp_pe(),'cost ',d3,d0,d7 write(0,'(i3,a,3f)') mpp_pe(),'cost ',d4,d5,d6 write(0,'(i3,a,2f)') mpp_pe(),'x,y ',x,y write(0,'(i3,a,4i5)') mpp_pe(),'io,jo ',io,jo,di,dj write(0,'(i3,a,2i5,a,2i3)') mpp_pe(),'diamonds, find_cell_by_search: false negative 2 i,j=',i-is,j-js,' di,dj=',di,dj write(0,'(i3,a,2i5,a,2f8.3)') mpp_pe(),'diamonds, find_cell_by_search: false negative 2 io,jo=',io,jo write(0,'(i3,a,2i5,a,2f8.3)') mpp_pe(),'diamonds, find_cell_by_search: false negative 2 i,j=',i,j,' targ=',x,y return endif find_cell_by_search=.false. contains ! # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # real function dcost(x1, y1, x2, y2) ! Arguments real, intent(in) :: x1, x2, y1, y2 ! Local variables real :: x1m x1m=modulo(x1-(x2-180.),360.)+(x2-180.) ! dcost=(x2-x1)**2+(y2-y1)**2 dcost=(x2-x1m)**2+(y2-y1)**2 end function dcost ! # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # logical function find_better_min(grd, x, y, w, oi, oj) ! Arguments type(icebergs_gridded), intent(in) :: grd real, intent(in) :: x, y integer, intent(in) :: w integer, intent(inout) :: oi, oj ! Local variables integer :: i,j,xs,xe,ys,ye real :: dmin, dcst xs=max(grd%isc, oi-w) xe=min(grd%iec, oi+w) ys=max(grd%jsc, oj-w) ye=min(grd%jec, oj+w) find_better_min=.false. dmin=dcost(x,y,grd%lonc(oi,oj),grd%latc(oi,oj)) do j=ys,ye; do i=xs,xe dcst=dcost(x,y,grd%lonc(i,j),grd%latc(i,j)) if (dcstgrd%ied) return if (j>grd%jed) return x1=grd%lon(i-1,j-1) y1=grd%lat(i-1,j-1) x2=grd%lon(i ,j-1) y2=grd%lat(i ,j-1) x3=grd%lon(i ,j ) y3=grd%lat(i ,j ) x4=grd%lon(i-1,j ) y4=grd%lat(i-1,j ) if (present(explain).and.explain) then write(stderrunit,'(a,4f12.6)') 'pos_within_cell: x1..x4 ',x1,x2,x3,x4 write(stderrunit,'(a,2f12.6)') 'pos_within_cell: x',x write(stderrunit,'(a,4f12.6)') 'pos_within_cell: y1..y4 ',y1,y2,y3,y4 write(stderrunit,'(a,2f12.6)') 'pos_within_cell: y',y endif if (max(y1,y2,y3,y4)<89.999) then call calc_xiyj(x1, x2, x3, x4, y1, y2, y3, y4, x, y, xi, yj, explain=explain) else if (debug) write(stderrunit,*) 'diamonds, pos_within_cell: working in tangential plane!' xx=(90.-y)*cos(x*pi_180) yy=(90.-y)*sin(x*pi_180) x1=(90.-y1)*cos(grd%lon(i-1,j-1)*pi_180) y1=(90.-y1)*sin(grd%lon(i-1,j-1)*pi_180) x2=(90.-y2)*cos(grd%lon(i ,j-1)*pi_180) y2=(90.-y2)*sin(grd%lon(i ,j-1)*pi_180) x3=(90.-y3)*cos(grd%lon(i ,j )*pi_180) y3=(90.-y3)*sin(grd%lon(i ,j )*pi_180) x4=(90.-y4)*cos(grd%lon(i-1,j )*pi_180) y4=(90.-y4)*sin(grd%lon(i-1,j )*pi_180) if (present(explain).and.explain) then write(stderrunit,'(a,4f12.6)') 'pos_within_cell: x1..x4 ',x1,x2,x3,x4 write(stderrunit,'(a,2f12.6)') 'pos_within_cell: x',xx write(stderrunit,'(a,4f12.6)') 'pos_within_cell: y1..y4 ',y1,y2,y3,y4 write(stderrunit,'(a,2f12.6)') 'pos_within_cell: y',yy endif call calc_xiyj(x1, x2, x3, x4, y1, y2, y3, y4, xx, yy, xi, yj, explain=explain) if (is_point_in_cell(grd, x, y, i, j)) then if (abs(xi-0.5)>0.5.or.abs(yj-0.5)>0.5) then ! Scale internal coordinates to be consistent with is_point_in_cell() ! Note: this is intended to fix the inconsistency between the tangent plane ! and lat-lon calculations fac=2.*max( abs(xi-0.5), abs(yj-0.5) ); fac=max(1., fac) xi=0.5+(xi-0.5)/fac yj=0.5+(yj-0.5)/fac if (debug) call error_mesg('diamonds, pos_within_cell', 'in cell so scaling internal coordinates!', WARNING) endif else if (abs(xi-0.5)<=0.5.and.abs(yj-0.5)<=0.5) then if (debug) call error_mesg('diamonds, pos_within_cell', 'out of cell but coordinates <=0.5!', WARNING) endif endif endif if (present(explain).and.explain) then write(stderrunit,'(a,2f12.6)') 'pos_within_cell: xi,yj=',xi,yj endif !if (.not. is_point_in_cell(grd, x, y, i, j) ) then ! write(stderrunit,'(a,i3,a,8f8.2,a)') 'diamonds, pos_within_cell: (',mpp_pe(),') ', & ! x1, y1, x2, y2, x3, y3, x4, y4, ' NOT IN CELL!' !endif if (xi.ge.0. .and. xi.le.1. .and. yj.ge.0. .and. yj.le.1.) then pos_within_cell=is_point_in_cell(grd, x, y, i, j, explain=explain) if (.not. pos_within_cell .and. verbose) then if (debug) call error_mesg('diamonds, pos_within_cell', 'pos_within_cell is in cell BUT is_point_in_cell disagrees!', WARNING) endif !pos_within_cell=.true. ! commenting this out makes pos_within_cell agree with is_point_in_cell endif contains subroutine calc_xiyj(x1, x2, x3, x4, y1, y2, y3, y4, x, y, xi, yj, explain) ! Arguments real, intent(in) :: x1, x2, x3, x4, y1, y2, y3, y4, x, y real, intent(out) :: xi, yj logical, intent(in), optional :: explain ! Local variables real :: alpha, beta, gamma, delta, epsilon, kappa, a, b, c, d, dx, dy, yy1, yy2 logical :: expl=.false. integer :: stderrunit stderrunit=stderr() expl=.false. if (present(explain).and.explain) expl=.true. alpha=x2-x1 delta=y2-y1 beta=x4-x1 epsilon=y4-y1 gamma=(x3-x1)-(alpha+beta) kappa=(y3-y1)-(delta+epsilon) if (expl) write(stderrunit,'(a,1p6e12.4)') 'calc_xiyj: coeffs alpha,beta,gamma',alpha,beta,gamma,delta,epsilon,kappa if (expl) write(stderrunit,'(a,1p6e12.4)') 'calc_xiyj: coeffs delta,epsilon,kappa',alpha,beta,gamma,delta,epsilon,kappa a=(kappa*beta-gamma*epsilon) dx=modulo(x-(x1-180.),360.)+(x1-180.)-x1 dy=y-y1 b=(delta*beta-alpha*epsilon)-(kappa*dx-gamma*dy) c=(alpha*dy-delta*dx) if (expl) write(stderrunit,'(a,1p6e12.4)') 'calc_xiyj: coeffs dx,dy=',dx,dy if (expl) write(stderrunit,'(a,1p6e12.4)') 'calc_xiyj: coeffs A,B,C=',a,b,c if (abs(a)>1.e-12) then d=0.25*(b**2)-a*c if (expl) write(stderrunit,'(a,1p6e12.4)') 'calc_xiyj: coeffs D=',d if (d.ge.0.) then if (expl) write(stderrunit,'(a,1p3e12.4)') 'Roots for b/2a, sqrt(d) = ',-0.5*b/a,sqrt(d)/a yy1=-(0.5*b+sqrt(d))/a yy2=-(0.5*b-sqrt(d))/a if (abs(yy1-0.5).lt.abs(yy2-0.5)) then; yj=yy1; else; yj=yy2; endif if (expl) write(stderrunit,'(a,1p3e12.4)') 'Roots for y = ',yy1,yy2,yj else write(stderrunit,'(a,i3,4f8.2)') 'calc_xiyj: x1..x4 ',mpp_pe(),x1,x2,x3,x4 write(stderrunit,'(a,i3,3f8.2)') 'calc_xiyj: x2..x4 - x1',mpp_pe(),x2-x1,x3-x1,x4-x1 write(stderrunit,'(a,i3,4f8.2)') 'calc_xiyj: y1..y4 ',mpp_pe(),y1,y2,y3,y4 write(stderrunit,'(a,i3,3f8.2)') 'calc_xiyj: y2..y4 - x1',mpp_pe(),y2-y1,y3-y1,y4-y1 write(stderrunit,'(a,i3,1p6e12.4)') 'calc_xiyj: coeffs alpha..kappa',mpp_pe(),alpha,beta,gamma,delta,epsilon,kappa write(stderrunit,'(a,i3)') 'calc_xiyj: b<0 in quadratic root solver!!!!',mpp_pe() write(stderrunit,'(a,i3,1p6e12.4)') 'calc_xiyj: coeffs a,b,c,d,dx,dy',mpp_pe(),a,b,c,d,dx,dy call error_mesg('diamonds, calc_xiyj', 'We have complex roots. The grid must be very distorted!', FATAL) endif else if (b.ne.0.) then yj=-c/b else yj=0. endif endif a=(alpha+gamma*yj) b=(delta+kappa*yj) if (a.ne.0.) then xi=(dx-beta*yj)/a elseif (b.ne.0.) then xi=(dy-epsilon*yj)/b else c=(epsilon*alpha-beta*delta)+(epsilon*gamma-beta*kappa)*yj if (c.ne.0.) then xi=(epsilon*dx-beta*dy)/c else write(stderrunit,'(a,i3,4f8.2)') 'calc_xiyj: x1..x4 ',mpp_pe(),x1,x2,x3,x4 write(stderrunit,'(a,i3,3f8.2)') 'calc_xiyj: x2..x4 - x1',mpp_pe(),x2-x1,x3-x1,x4-x1 write(stderrunit,'(a,i3,4f8.2)') 'calc_xiyj: y1..y4 ',mpp_pe(),y1,y2,y3,y4 write(stderrunit,'(a,i3,3f8.2)') 'calc_xiyj: y2..y4 - x1',mpp_pe(),y2-y1,y3-y1,y4-y1 write(stderrunit,'(a,i3,1p6e12.4)') 'calc_xiyj: coeffs alpha..kappa',mpp_pe(),alpha,beta,gamma,delta,epsilon,kappa write(stderrunit,'(a,i3,1p2e12.4)') 'calc_xiyj: coeffs a,b',mpp_pe(),a,b call error_mesg('diamonds, calc_xiyj', 'Can not invert either linear equaton for xi! This should not happen!', FATAL) endif endif if (expl) write(stderrunit,'(a,2e12.4)') 'calc_xiyj: xi,yj=',xi,yj end subroutine calc_xiyj end function pos_within_cell ! ############################################################################## subroutine check_position(grd, berg, label) ! Arguments type(icebergs_gridded), pointer :: grd type(iceberg), pointer :: berg character(len=*) :: label ! Local variables real :: xi, yj logical :: lret integer :: stderrunit stderrunit=stderr() lret=pos_within_cell(grd, berg%lon, berg%lat, berg%ine, berg%jne, xi, yj) if (xi.ne.berg%xi.or.yj.ne.berg%yj) then write(stderrunit,'("diamonds: check_position (",i4,") b%x,x,-=",3(es12.4,x),a)'),mpp_pe(),berg%xi,xi,berg%xi-xi,label write(stderrunit,'("diamonds: check_position (",i4,") b%y,y,-=",3(es12.4,x),a)'),mpp_pe(),berg%yj,yj,berg%yj-yj,label call print_berg(stderrunit, berg, 'check_position') call error_mesg('diamonds, check_position','berg has inconsistent xi,yj!',FATAL) endif end subroutine check_position ! ############################################################################## real function sum_mass(first,justbits,justbergs) ! Arguments type(iceberg), pointer :: first logical, intent(in), optional :: justbits, justbergs ! Local variables type(iceberg), pointer :: this sum_mass=0. this=>first do while(associated(this)) if (present(justbergs)) then sum_mass=sum_mass+this%mass*this%mass_scaling elseif (present(justbits)) then sum_mass=sum_mass+this%mass_of_bits*this%mass_scaling else sum_mass=sum_mass+(this%mass+this%mass_of_bits)*this%mass_scaling endif this=>this%next enddo end function sum_mass ! ############################################################################## real function sum_heat(first,justbits,justbergs) ! Arguments type(iceberg), pointer :: first logical, intent(in), optional :: justbits, justbergs ! Local variables type(iceberg), pointer :: this real :: dm sum_heat=0. this=>first do while(associated(this)) dm=0. if (present(justbergs)) then dm=this%mass*this%mass_scaling elseif (present(justbits)) then dm=this%mass_of_bits*this%mass_scaling else dm=(this%mass+this%mass_of_bits)*this%mass_scaling endif sum_heat=sum_heat+dm*this%heat_density this=>this%next enddo end function sum_heat ! ############################################################################## subroutine icebergs_stock_pe(bergs, index, value) ! Modules use stock_constants_mod, only : ISTOCK_WATER, ISTOCK_HEAT ! Arguments type(icebergs), pointer :: bergs integer, intent(in) :: index real, intent(out) :: value ! Local variables type(icebergs_gridded), pointer :: grd real :: berg_mass, stored_mass ! For convenience grd=>bergs%grd select case (index) case (ISTOCK_WATER) berg_mass=sum_mass(bergs%first) stored_mass=sum( grd%stored_ice(grd%isc:grd%iec,grd%jsc:grd%jec,:) ) value=stored_mass+berg_mass case (ISTOCK_HEAT) berg_mass=sum_mass(bergs%first) stored_mass=sum( grd%stored_ice(grd%isc:grd%iec,grd%jsc:grd%jec,:) ) value=-(stored_mass+berg_mass)*HLF ! HLF is in (J/kg) from constants_mod case default value = 0.0 end select end subroutine icebergs_stock_pe ! ############################################################################## subroutine icebergs_save_restart(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables if (.not.associated(bergs)) return call mpp_clock_begin(bergs%clock_iow) call bergs_chksum(bergs, 'write_restart bergs') call write_restart(bergs) call mpp_clock_end(bergs%clock_iow) end subroutine icebergs_save_restart ! ############################################################################## subroutine icebergs_end(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables type(iceberg), pointer :: this, next if (.not.associated(bergs)) return call icebergs_save_restart(bergs) call mpp_clock_begin(bergs%clock_ini) ! Delete bergs and structures this=>bergs%first do while (associated(this)) next=>this%next call move_trajectory(bergs, this) call destroy_iceberg(this) this=>next enddo if (associated(bergs%trajectories)) then call write_trajectory(bergs%trajectories) endif deallocate(bergs%grd%lon) deallocate(bergs%grd%lat) deallocate(bergs%grd%lonc) deallocate(bergs%grd%latc) deallocate(bergs%grd%dx) deallocate(bergs%grd%dy) deallocate(bergs%grd%area) deallocate(bergs%grd%msk) deallocate(bergs%grd%cos) deallocate(bergs%grd%sin) deallocate(bergs%grd%calving) deallocate(bergs%grd%calving_hflx) deallocate(bergs%grd%stored_heat) deallocate(bergs%grd%floating_melt) deallocate(bergs%grd%berg_melt) deallocate(bergs%grd%melt_buoy) deallocate(bergs%grd%melt_eros) deallocate(bergs%grd%melt_conv) deallocate(bergs%grd%bergy_src) deallocate(bergs%grd%bergy_melt) deallocate(bergs%grd%bergy_mass) deallocate(bergs%grd%virtual_area) deallocate(bergs%grd%mass) deallocate(bergs%grd%mass_on_ocean) deallocate(bergs%grd%tmp) deallocate(bergs%grd%tmpc) deallocate(bergs%grd%stored_ice) deallocate(bergs%grd%real_calving) deallocate(bergs%grd%uo) deallocate(bergs%grd%vo) deallocate(bergs%grd%ui) deallocate(bergs%grd%vi) deallocate(bergs%grd%ua) deallocate(bergs%grd%va) deallocate(bergs%grd%ssh) deallocate(bergs%grd%sst) deallocate(bergs%grd%cn) deallocate(bergs%grd%hi) deallocate(bergs%grd%domain) deallocate(bergs%grd) deallocate(bergs%initial_mass) deallocate(bergs%distribution) deallocate(bergs%initial_thickness) deallocate(bergs%initial_width) deallocate(bergs%initial_length) call dealloc_buffer(bergs%obuffer_n) call dealloc_buffer(bergs%obuffer_s) call dealloc_buffer(bergs%obuffer_e) call dealloc_buffer(bergs%obuffer_w) call dealloc_buffer(bergs%ibuffer_n) call dealloc_buffer(bergs%ibuffer_s) call dealloc_buffer(bergs%ibuffer_e) call dealloc_buffer(bergs%ibuffer_w) call mpp_clock_end(bergs%clock_ini) deallocate(bergs) if (mpp_pe()==mpp_root_pe()) write(*,'(a,i)') 'diamonds: icebergs_end complete',mpp_pe() contains subroutine dealloc_buffer(buff) ! Arguments type(buffer), pointer :: buff ! Local variables if (associated(buff)) then if (associated(buff%data)) deallocate(buff%data) deallocate(buff) endif end subroutine dealloc_buffer end subroutine icebergs_end ! ############################################################################## subroutine invert_tau_for_du(u,v) ! Arguments real, dimension(:,:),intent(inout) :: u, v ! Local variables integer :: i, j real :: cd, cddvmod, tau2 cd=0.0015 do j=lbound(u,2), ubound(u,2) do i=lbound(u,1), ubound(u,1) tau2=u(i,j)*u(i,j)+v(i,j)*v(i,j) cddvmod=sqrt(cd*sqrt(tau2)) if (cddvmod.ne.0.) then u(i,j)=u(i,j)/cddvmod v(i,j)=v(i,j)/cddvmod else u(i,j)=0. v(i,j)=0. endif enddo enddo end subroutine invert_tau_for_du ! ############################################################################## subroutine sanitize_field(arr,val) ! Arguments real, dimension(:,:),intent(inout) :: arr real, intent(in) :: val ! Local variables integer :: i, j do j=lbound(arr,2), ubound(arr,2) do i=lbound(arr,1), ubound(arr,1) if (abs(arr(i,j)).ge.val) arr(i,j)=0. enddo enddo end subroutine sanitize_field ! ############################################################################## subroutine write_restart(bergs) ! Arguments type(icebergs), pointer :: bergs ! Local variables integer :: iret, ncid, i_dim, i integer :: lonid, latid, uvelid, vvelid, ineid, jneid integer :: massid, thicknessid, lengthid, widthid integer :: start_lonid, start_latid, start_yearid, start_dayid, start_massid integer :: scaling_id, mass_of_bits_id, heat_density_id character(len=28) :: filename type(iceberg), pointer :: this integer :: stderrunit stderrunit=stderr() ! Only create a restart file for this PE if we have anything to say if (associated(bergs%first)) then write(filename(1:28),'("RESTART/icebergs.res.nc.",I4.4)') mpp_pe() if (verbose) write(*,'(2a)') 'diamonds, write_restart: creating ',filename iret = nf_create(filename, NF_CLOBBER, ncid) if (iret .ne. NF_NOERR) write(stderrunit,*) 'diamonds, write_restart: nf_create failed' ! Dimensions iret = nf_def_dim(ncid, 'i', NF_UNLIMITED, i_dim) if (iret .ne. NF_NOERR) write(stderrunit,*) 'diamonds, write_restart: nf_def_dim i failed' ! Variables lonid = def_var(ncid, 'lon', NF_DOUBLE, i_dim) latid = def_var(ncid, 'lat', NF_DOUBLE, i_dim) uvelid = def_var(ncid, 'uvel', NF_DOUBLE, i_dim) vvelid = def_var(ncid, 'vvel', NF_DOUBLE, i_dim) massid = def_var(ncid, 'mass', NF_DOUBLE, i_dim) ineid = def_var(ncid, 'ine', NF_INT, i_dim) jneid = def_var(ncid, 'jne', NF_INT, i_dim) thicknessid = def_var(ncid, 'thickness', NF_DOUBLE, i_dim) widthid = def_var(ncid, 'width', NF_DOUBLE, i_dim) lengthid = def_var(ncid, 'length', NF_DOUBLE, i_dim) start_lonid = def_var(ncid, 'start_lon', NF_DOUBLE, i_dim) start_latid = def_var(ncid, 'start_lat', NF_DOUBLE, i_dim) start_yearid = def_var(ncid, 'start_year', NF_INT, i_dim) start_dayid = def_var(ncid, 'start_day', NF_DOUBLE, i_dim) start_massid = def_var(ncid, 'start_mass', NF_DOUBLE, i_dim) scaling_id = def_var(ncid, 'mass_scaling', NF_DOUBLE, i_dim) mass_of_bits_id = def_var(ncid, 'mass_of_bits', NF_DOUBLE, i_dim) heat_density_id = def_var(ncid, 'heat_density', NF_DOUBLE, i_dim) ! Attributes call put_att(ncid, lonid, 'long_name', 'longitude') call put_att(ncid, lonid, 'units', 'degrees_E') call put_att(ncid, latid, 'long_name', 'latitude') call put_att(ncid, latid, 'units', 'degrees_N') call put_att(ncid, uvelid, 'long_name', 'zonal velocity') call put_att(ncid, uvelid, 'units', 'm/s') call put_att(ncid, vvelid, 'long_name', 'meridional velocity') call put_att(ncid, vvelid, 'units', 'm/s') call put_att(ncid, ineid, 'long_name', 'i index') call put_att(ncid, ineid, 'units', 'none') call put_att(ncid, jneid, 'long_name', 'j index') call put_att(ncid, jneid, 'units', 'none') call put_att(ncid, massid, 'long_name', 'mass') call put_att(ncid, massid, 'units', 'kg') call put_att(ncid, thicknessid, 'long_name', 'thickness') call put_att(ncid, thicknessid, 'units', 'm') call put_att(ncid, widthid, 'long_name', 'width') call put_att(ncid, widthid, 'units', 'm') call put_att(ncid, lengthid, 'long_name', 'length') call put_att(ncid, lengthid, 'units', 'm') call put_att(ncid, start_lonid, 'long_name', 'longitude of calving location') call put_att(ncid, start_lonid, 'units', 'degrees_E') call put_att(ncid, start_latid, 'long_name', 'latitude of calving location') call put_att(ncid, start_latid, 'units', 'degrees_N') call put_att(ncid, start_yearid, 'long_name', 'calendar year of calving event') call put_att(ncid, start_yearid, 'units', 'years') call put_att(ncid, start_dayid, 'long_name', 'year day of calving event') call put_att(ncid, start_dayid, 'units', 'days') call put_att(ncid, start_massid, 'long_name', 'initial mass of calving berg') call put_att(ncid, start_massid, 'units', 'kg') call put_att(ncid, scaling_id, 'long_name', 'scaling factor for mass of calving berg') call put_att(ncid, scaling_id, 'units', 'none') call put_att(ncid, mass_of_bits_id, 'long_name', 'mass of bergy bits') call put_att(ncid, mass_of_bits_id, 'units', 'kg') call put_att(ncid, heat_density_id, 'long_name', 'heat density') call put_att(ncid, heat_density_id, 'units', 'J/kg') iret = nf_put_att_int(ncid, NCGLOBAL, 'file_format_major_version', NF_INT, 1, file_format_major_version) iret = nf_put_att_int(ncid, NCGLOBAL, 'file_format_minor_version', NF_INT, 3, file_format_minor_version) ! End define mode iret = nf_enddef(ncid) ! Write variables !this=>bergs%first; i=0 this=>last_berg(bergs%first); i=0 do while (associated(this)) i=i+1 call put_double(ncid, lonid, i, this%lon) call put_double(ncid, latid, i, this%lat) call put_double(ncid, uvelid, i, this%uvel) call put_double(ncid, vvelid, i, this%vvel) call put_int(ncid, ineid, i, this%ine) call put_int(ncid, jneid, i, this%jne) call put_double(ncid, massid, i, this%mass) call put_double(ncid, thicknessid, i, this%thickness) call put_double(ncid, widthid, i, this%width) call put_double(ncid, lengthid, i, this%length) call put_double(ncid, start_lonid, i, this%start_lon) call put_double(ncid, start_latid, i, this%start_lat) call put_int(ncid, start_yearid, i, this%start_year) call put_double(ncid, start_dayid, i, this%start_day) call put_double(ncid, start_massid, i, this%start_mass) call put_double(ncid, scaling_id, i, this%mass_scaling) call put_double(ncid, mass_of_bits_id, i, this%mass_of_bits) call put_double(ncid, heat_density_id, i, this%heat_density) !this=>this%next this=>this%prev enddo ! Finish up iret = nf_close(ncid) if (iret .ne. NF_NOERR) write(stderrunit,*) 'diamonds, write_restart: nf_close failed' endif ! associated(bergs%first) ! Write stored ice filename='RESTART/calving.res.nc' if (verbose.and.mpp_pe().eq.mpp_root_pe()) write(stderrunit,'(2a)') 'diamonds, write_restart: writing ',filename call grd_chksum3(bergs%grd, bergs%grd%stored_ice, 'write stored_ice') call write_data(filename, 'stored_ice', bergs%grd%stored_ice, bergs%grd%domain) call grd_chksum2(bergs%grd, bergs%grd%stored_heat, 'write stored_heat') call write_data(filename, 'stored_heat', bergs%grd%stored_heat, bergs%grd%domain) contains function last_berg(berg) ! Arguments type(iceberg), pointer :: last_berg, berg ! Local variables last_berg=>berg do while (associated(last_berg%next)) last_berg=>last_berg%next enddo end function last_berg end subroutine write_restart ! ############################################################################## subroutine write_trajectory(trajectory) ! Arguments type(xyt), pointer :: trajectory ! Local variables integer :: iret, ncid, i_dim, i integer :: lonid, latid, yearid, dayid, uvelid, vvelid integer :: uoid, void, uiid, viid, uaid, vaid, sshxid, sshyid, sstid integer :: cnid, hiid integer :: mid, did, wid, lid, mbid, hdid character(len=30) :: filename type(xyt), pointer :: this, next integer :: stderrunit stderrunit=stderr() write(filename(1:30),'("iceberg_trajectories.nc.",I4.4)') mpp_pe() if (debug) write(stderrunit,*) 'diamonds, write_trajectory: creating ',filename iret = nf_create(filename, NF_CLOBBER, ncid) if (iret .ne. NF_NOERR) write(stderrunit,*) 'diamonds, write_trajectory: nf_create failed' ! Dimensions iret = nf_def_dim(ncid, 'i', NF_UNLIMITED, i_dim) if (iret .ne. NF_NOERR) write(stderrunit,*) 'diamonds, write_trajectory: nf_def_dim i failed' ! Variables lonid = def_var(ncid, 'lon', NF_DOUBLE, i_dim) latid = def_var(ncid, 'lat', NF_DOUBLE, i_dim) yearid = def_var(ncid, 'year', NF_INT, i_dim) dayid = def_var(ncid, 'day', NF_DOUBLE, i_dim) uvelid = def_var(ncid, 'uvel', NF_DOUBLE, i_dim) vvelid = def_var(ncid, 'vvel', NF_DOUBLE, i_dim) uoid = def_var(ncid, 'uo', NF_DOUBLE, i_dim) void = def_var(ncid, 'vo', NF_DOUBLE, i_dim) uiid = def_var(ncid, 'ui', NF_DOUBLE, i_dim) viid = def_var(ncid, 'vi', NF_DOUBLE, i_dim) uaid = def_var(ncid, 'ua', NF_DOUBLE, i_dim) vaid = def_var(ncid, 'va', NF_DOUBLE, i_dim) mid = def_var(ncid, 'mass', NF_DOUBLE, i_dim) mbid = def_var(ncid, 'mass_of_bits', NF_DOUBLE, i_dim) hdid = def_var(ncid, 'heat_density', NF_DOUBLE, i_dim) did = def_var(ncid, 'thickness', NF_DOUBLE, i_dim) wid = def_var(ncid, 'width', NF_DOUBLE, i_dim) lid = def_var(ncid, 'length', NF_DOUBLE, i_dim) sshxid = def_var(ncid, 'ssh_x', NF_DOUBLE, i_dim) sshyid = def_var(ncid, 'ssh_y', NF_DOUBLE, i_dim) sstid = def_var(ncid, 'sst', NF_DOUBLE, i_dim) cnid = def_var(ncid, 'cn', NF_DOUBLE, i_dim) hiid = def_var(ncid, 'hi', NF_DOUBLE, i_dim) ! Attributes iret = nf_put_att_int(ncid, NCGLOBAL, 'file_format_major_version', NF_INT, 1, 0) iret = nf_put_att_int(ncid, NCGLOBAL, 'file_format_minor_version', NF_INT, 1, 1) call put_att(ncid, lonid, 'long_name', 'longitude') call put_att(ncid, lonid, 'units', 'degrees_E') call put_att(ncid, latid, 'long_name', 'latitude') call put_att(ncid, latid, 'units', 'degrees_N') call put_att(ncid, yearid, 'long_name', 'year') call put_att(ncid, yearid, 'units', 'years') call put_att(ncid, dayid, 'long_name', 'year day') call put_att(ncid, dayid, 'units', 'days') call put_att(ncid, uvelid, 'long_name', 'zonal spped') call put_att(ncid, uvelid, 'units', 'm/s') call put_att(ncid, vvelid, 'long_name', 'meridional spped') call put_att(ncid, vvelid, 'units', 'm/s') call put_att(ncid, uoid, 'long_name', 'ocean zonal spped') call put_att(ncid, uoid, 'units', 'm/s') call put_att(ncid, void, 'long_name', 'ocean meridional spped') call put_att(ncid, void, 'units', 'm/s') call put_att(ncid, uiid, 'long_name', 'ice zonal spped') call put_att(ncid, uiid, 'units', 'm/s') call put_att(ncid, viid, 'long_name', 'ice meridional spped') call put_att(ncid, viid, 'units', 'm/s') call put_att(ncid, uaid, 'long_name', 'atmos zonal spped') call put_att(ncid, uaid, 'units', 'm/s') call put_att(ncid, vaid, 'long_name', 'atmos meridional spped') call put_att(ncid, vaid, 'units', 'm/s') call put_att(ncid, mid, 'long_name', 'mass') call put_att(ncid, mid, 'units', 'kg') call put_att(ncid, mbid, 'long_name', 'mass_of_bits') call put_att(ncid, mbid, 'units', 'kg') call put_att(ncid, hdid, 'long_name', 'heat_density') call put_att(ncid, hdid, 'units', 'J/kg') call put_att(ncid, did, 'long_name', 'thickness') call put_att(ncid, did, 'units', 'm') call put_att(ncid, wid, 'long_name', 'width') call put_att(ncid, wid, 'units', 'm') call put_att(ncid, lid, 'long_name', 'length') call put_att(ncid, lid, 'units', 'm') call put_att(ncid, sshxid, 'long_name', 'sea surface height gradient_x') call put_att(ncid, sshxid, 'units', 'non-dim') call put_att(ncid, sshyid, 'long_name', 'sea surface height gradient_y') call put_att(ncid, sshyid, 'units', 'non-dim') call put_att(ncid, sstid, 'long_name', 'sea surface temperature') call put_att(ncid, sstid, 'units', 'degrees_C') call put_att(ncid, cnid, 'long_name', 'sea ice concentration') call put_att(ncid, cnid, 'units', 'none') call put_att(ncid, hiid, 'long_name', 'sea ice thickness') call put_att(ncid, hiid, 'units', 'm') ! End define mode iret = nf_enddef(ncid) ! Write variables this=>trajectory; i=0 do while (associated(this)) i=i+1 call put_double(ncid, lonid, i, this%lon) call put_double(ncid, latid, i, this%lat) call put_int(ncid, yearid, i, this%year) call put_double(ncid, dayid, i, this%day) call put_double(ncid, uvelid, i, this%uvel) call put_double(ncid, vvelid, i, this%vvel) call put_double(ncid, uoid, i, this%uo) call put_double(ncid, void, i, this%vo) call put_double(ncid, uiid, i, this%ui) call put_double(ncid, viid, i, this%vi) call put_double(ncid, uaid, i, this%ua) call put_double(ncid, vaid, i, this%va) call put_double(ncid, mid, i, this%mass) call put_double(ncid, hdid, i, this%heat_density) call put_double(ncid, did, i, this%thickness) call put_double(ncid, wid, i, this%width) call put_double(ncid, lid, i, this%length) call put_double(ncid, sshxid, i, this%ssh_x) call put_double(ncid, sshyid, i, this%ssh_y) call put_double(ncid, sstid, i, this%sst) call put_double(ncid, cnid, i, this%cn) call put_double(ncid, hiid, i, this%hi) next=>this%next deallocate(this) this=>next enddo trajectory=>NULL() ! Finish up iret = nf_close(ncid) if (iret .ne. NF_NOERR) write(stderrunit,*) 'diamonds, write_trajectory: nf_close failed',mpp_pe(),filename end subroutine write_trajectory ! ############################################################################## integer function inq_var(ncid, var, unsafe) ! Arguments integer, intent(in) :: ncid character(len=*), intent(in) :: var logical, optional, intent(in) :: unsafe ! Local variables integer :: iret integer :: stderrunit stderrunit=stderr() iret=nf_inq_varid(ncid, var, inq_var) if (iret .ne. NF_NOERR) then if (.not.present(unsafe) .or. (present(unsafe).and..not.unsafe)) then write(stderrunit,*) 'diamonds, inq_var: nf_inq_varid ',var,' failed' call error_mesg('diamonds, inq_var', 'netcdf function returned a failure!', FATAL) else inq_var=-1 endif endif end function inq_var ! ############################################################################## integer function def_var(ncid, var, ntype, idim) ! Arguments integer, intent(in) :: ncid, ntype, idim character(len=*), intent(in) :: var ! Local variables integer :: iret integer :: stderrunit stderrunit=stderr() iret = nf_def_var(ncid, var, ntype, 1, idim, def_var) if (iret .ne. NF_NOERR) then write(stderrunit,*) 'diamonds, def_var: nf_def_var failed for ',trim(var) call error_mesg('diamonds, def_var', 'netcdf function returned a failure!', FATAL) endif end function def_var ! ############################################################################## subroutine put_att(ncid, id, att, attval) ! Arguments integer, intent(in) :: ncid, id character(len=*), intent(in) :: att, attval ! Local variables integer :: vallen, iret integer :: stderrunit stderrunit=stderr() vallen=LEN_TRIM(attval) iret = nf_put_att_text(ncid, id, att, vallen, attval) if (iret .ne. NF_NOERR) then write(stderrunit,*) 'diamonds, put_att: nf_put_att_text failed adding', & trim(att),' = ',trim(attval) call error_mesg('diamonds, put_att', 'netcdf function returned a failure!', FATAL) endif end subroutine put_att ! ############################################################################## real function get_double(ncid, id, i) ! Arguments integer, intent(in) :: ncid, id, i ! Local variables integer :: iret integer :: stderrunit stderrunit=stderr() iret=nf_get_var1_double(ncid, id, i, get_double) if (iret .ne. NF_NOERR) then write(stderrunit,*) 'diamonds, get_double: nf_get_var1_double failed reading' call error_mesg('diamonds, get_double', 'netcdf function returned a failure!', FATAL) endif end function get_double ! ############################################################################## integer function get_int(ncid, id, i) ! Arguments integer, intent(in) :: ncid, id, i ! Local variables integer :: iret integer :: stderrunit stderrunit=stderr() iret=nf_get_var1_int(ncid, id, i, get_int) if (iret .ne. NF_NOERR) then write(stderrunit,*) 'diamonds, get_int: nf_get_var1_int failed reading' call error_mesg('diamonds, get_int', 'netcdf function returned a failure!', FATAL) endif end function get_int ! ############################################################################## subroutine put_double(ncid, id, i, val) ! Arguments integer, intent(in) :: ncid, id, i real, intent(in) :: val ! Local variables integer :: iret integer :: stderrunit stderrunit=stderr() iret = nf_put_vara_double(ncid, id, i, 1, val) if (iret .ne. NF_NOERR) then write(stderrunit,*) 'diamonds, put_double: nf_put_vara_double failed writing' call error_mesg('diamonds, put_double', 'netcdf function returned a failure!', FATAL) endif end subroutine put_double ! ############################################################################## subroutine put_int(ncid, id, i, val) ! Arguments integer, intent(in) :: ncid, id, i, val ! Local variables integer :: iret integer :: stderrunit stderrunit=stderr() iret = nf_put_vara_int(ncid, id, i, 1, val) if (iret .ne. NF_NOERR) then write(stderrunit,*) 'diamonds, put_int: nf_put_vara_int failed writing' call error_mesg('diamonds, put_int', 'netcdf function returned a failure!', FATAL) endif end subroutine put_int ! ############################################################################## subroutine print_fld(grd, fld, label) ! Arguments type(icebergs_gridded), pointer :: grd real, intent(in) :: fld(grd%isd:grd%ied,grd%jsd:grd%jed) character(len=*) :: label ! Local variables integer :: i, j integer :: stderrunit stderrunit=stderr() write(stderrunit,'("pe=",i3,x,a8,32i10)') mpp_pe(),label,(i,i=grd%isd,grd%ied) do j=grd%jed,grd%jsd,-1 write(stderrunit,'("pe=",i3,x,i8,32es10.2)') mpp_pe(),j,(fld(i,j),i=grd%isd,grd%ied) enddo end subroutine print_fld ! ############################################################################## subroutine checksum_gridded(grd, label) ! Arguments type(icebergs_gridded), pointer :: grd character(len=*) :: label ! Local variables if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') 'diamonds: checksumming gridded data @ ',trim(label) ! external forcing call grd_chksum2(grd, grd%uo, 'uo') call grd_chksum2(grd, grd%vo, 'vo') call grd_chksum2(grd, grd%ua, 'ua') call grd_chksum2(grd, grd%va, 'va') call grd_chksum2(grd, grd%ui, 'ui') call grd_chksum2(grd, grd%vi, 'vi') call grd_chksum2(grd, grd%ssh, 'ssh') call grd_chksum2(grd, grd%sst, 'sst') call grd_chksum2(grd, grd%hi, 'hi') call grd_chksum2(grd, grd%cn, 'cn') call grd_chksum2(grd, grd%calving, 'calving') call grd_chksum2(grd, grd%calving_hflx, 'calving_hflx') ! state call grd_chksum2(grd, grd%mass, 'mass') call grd_chksum3(grd, grd%mass_on_ocean, 'mass_on_ocean') call grd_chksum3(grd, grd%stored_ice, 'stored_ice') call grd_chksum2(grd, grd%stored_heat, 'stored_heat') call grd_chksum2(grd, grd%melt_buoy, 'melt_b') call grd_chksum2(grd, grd%melt_eros, 'melt_e') call grd_chksum2(grd, grd%melt_conv, 'melt_v') call grd_chksum2(grd, grd%bergy_src, 'bergy_src') call grd_chksum2(grd, grd%bergy_melt, 'bergy_melt') call grd_chksum2(grd, grd%bergy_mass, 'bergy_mass') call grd_chksum2(grd, grd%virtual_area, 'varea') call grd_chksum2(grd, grd%floating_melt, 'floating_melt') call grd_chksum2(grd, grd%berg_melt, 'berg_melt') ! static call grd_chksum2(grd, grd%lon, 'lon') call grd_chksum2(grd, grd%lat, 'lat') call grd_chksum2(grd, grd%lonc, 'lonc') call grd_chksum2(grd, grd%latc, 'latc') call grd_chksum2(grd, grd%dx, 'dx') call grd_chksum2(grd, grd%dy, 'dy') call grd_chksum2(grd, grd%area, 'area') call grd_chksum2(grd, grd%msk, 'msk') call grd_chksum2(grd, grd%cos, 'cos') call grd_chksum2(grd, grd%sin, 'sin') end subroutine checksum_gridded ! ############################################################################## subroutine grd_chksum3(grd, fld, txt) ! Arguments type(icebergs_gridded), pointer :: grd real, dimension(:,:,:), intent(in) :: fld character(len=*), intent(in) :: txt ! Local variables integer :: i, j, k, halo, icount real :: mean, rms, SD, minv, maxv real, dimension(lbound(fld,1):ubound(fld,1), lbound(fld,2):ubound(fld,2), lbound(fld,3):ubound(fld,3)) :: tmp halo=grd%halo mean=0. rms=0. sd=0. icount=0 i=lbound(fld,1)+halo j=lbound(fld,2)+halo k=lbound(fld,3) minv=fld(i,j,k) maxv=fld(i,j,k) tmp(:,:,:)=0. do k=lbound(fld,3), ubound(fld,3) do j=lbound(fld,2)+halo, ubound(fld,2)-halo do i=lbound(fld,1)+halo, ubound(fld,1)-halo icount=icount+1 mean=mean+fld(i,j,k) rms=rms+fld(i,j,k)**2 minv=min(minv,fld(i,j,k)) maxv=max(maxv,fld(i,j,k)) tmp(i,j,k)=fld(i,j,k)*float(i+2*j+3*(k-1)) enddo enddo enddo call mpp_sum(icount) call mpp_sum(mean) call mpp_sum(rms) call mpp_min(minv) call mpp_max(maxv) mean=mean/float(icount) rms=sqrt(rms/float(icount)) do k=lbound(fld,3), ubound(fld,3) do j=lbound(fld,2)+halo, ubound(fld,2)-halo do i=lbound(fld,1)+halo, ubound(fld,1)-halo sd=sd+(fld(i,j,k)-mean)**2 enddo enddo enddo call mpp_sum(sd) sd=sqrt(sd/float(icount)) i=mpp_chksum( fld(lbound(fld,1)+halo:ubound(fld,1)-halo, & lbound(fld,2)+halo:ubound(fld,2)-halo,:) ) j=mpp_chksum( tmp(lbound(fld,1)+halo:ubound(fld,1)-halo, & lbound(fld,2)+halo:ubound(fld,2)-halo,:) ) if (mpp_pe().eq.mpp_root_pe()) & write(*,'("diamonds, grd_chksum3: ",a18,2(x,a,"=",i),5(x,a,"=",es16.9))') & txt, 'chksum', i, 'chksum2', j, 'min', minv, 'max', maxv, 'mean', mean, 'rms', rms, 'sd', sd #ifdef CHECKSUM_HALOS i=mpp_chksum( fld(lbound(fld,1):ubound(fld,1), & lbound(fld,2):ubound(fld,2),:) ) j=mpp_chksum( tmp(lbound(fld,1):ubound(fld,1), & lbound(fld,2):ubound(fld,2),:) ) if (mpp_pe().eq.mpp_root_pe()) & write(*,'("diamonds, grd_chksum3* ",a18,2(x,a,"=",i),5(x,a,"=",es16.9))') & txt, 'chksum', i, 'chksum2', j, 'min', minv, 'max', maxv, 'mean', mean, 'rms', rms, 'sd', sd #endif end subroutine grd_chksum3 ! ############################################################################## subroutine grd_chksum2(grd, fld, txt) ! Arguments type(icebergs_gridded), pointer :: grd real, dimension(grd%isd:grd%ied,grd%jsd:grd%jed), intent(in) :: fld character(len=*), intent(in) :: txt ! Local variables integer :: i, j, icount real :: mean, rms, SD, minv, maxv grd%tmp(:,:)=0. mean=0. rms=0. sd=0. icount=0 minv=fld(grd%isc,grd%jsc) maxv=fld(grd%isc,grd%jsc) do j=grd%jsd, grd%jed do i=grd%isd, grd%ied icount=icount+1 mean=mean+fld(i,j) rms=rms+fld(i,j)**2 minv=min(minv,fld(i,j)) maxv=max(maxv,fld(i,j)) grd%tmp(i,j)=fld(i,j)*float(i+2*j) enddo enddo call mpp_sum(icount) call mpp_sum(mean) call mpp_sum(rms) call mpp_min(minv) call mpp_max(maxv) mean=mean/float(icount) rms=sqrt(rms/float(icount)) do j=grd%jsc, grd%jec do i=grd%isc, grd%iec sd=sd+(fld(i,j)-mean)**2 enddo enddo call mpp_sum(sd) sd=sqrt(sd/float(icount)) i=mpp_chksum( fld(grd%isc:grd%iec,grd%jsc:grd%jec) ) j=mpp_chksum( grd%tmp(grd%isc:grd%iec,grd%jsc:grd%jec) ) if (mpp_pe().eq.mpp_root_pe()) & write(*,'("diamonds, grd_chksum2: ",a18,2(x,a,"=",i),5(x,a,"=",es16.9),x,a,"=",i)') & txt, 'chksum', i, 'chksum2', j, 'min', minv, 'max', maxv, 'mean', mean, 'rms', rms, 'sd', sd!, '#', icount #ifdef CHECKSUM_HALOS i=mpp_chksum( fld(grd%isd:grd%ied,grd%jsd:grd%jed) ) j=mpp_chksum( grd%tmp(grd%isd:grd%ied,grd%jsd:grd%jed) ) if (mpp_pe().eq.mpp_root_pe()) & write(*,'("diamonds, grd_chksum2* ",a18,2(x,a,"=",i),5(x,a,"=",es16.9),x,a,"=",i)') & txt, 'chksum', i, 'chksum2', j, 'min', minv, 'max', maxv, 'mean', mean, 'rms', rms, 'sd', sd!, '#', icount #endif end subroutine grd_chksum2 ! ############################################################################## subroutine bergs_chksum(bergs, txt, ignore_halo_violation) ! Arguments type(icebergs), pointer :: bergs character(len=*), intent(in) :: txt logical, optional :: ignore_halo_violation ! Local variables integer :: i, nbergs, ichk1, ichk2, ichk3, ichk4, ichk5, iberg real, allocatable :: fld(:,:), fld2(:,:) integer, allocatable :: icnt(:,:) type(iceberg), pointer :: this type(icebergs_gridded), pointer :: grd logical :: check_halo ! For convenience grd=>bergs%grd nbergs=count_bergs(bergs) call mpp_max(nbergs) allocate( fld( nbergs, 11 ) ) allocate( fld2( nbergs, 11 ) ) allocate( icnt( grd%isd:grd%ied, grd%jsd:grd%jed ) ) fld(:,:)=0. fld2(:,:)=0. icnt(:,:)=0 grd%tmp(:,:)=0. this=>bergs%first i=0; ichk4=0; ichk5=0 do while(associated(this)) i=i+1 iberg=berg_chksum(this) fld(i,1) = this%lon fld(i,2) = this%lat fld(i,3) = this%uvel fld(i,4) = this%vvel fld(i,5) = this%mass fld(i,6) = this%thickness fld(i,7) = this%width fld(i,8) = this%length fld(i,9) = time_hash(this) fld(i,10) = pos_hash(this) fld(i,11) = float(iberg) icnt(this%ine,this%jne)=icnt(this%ine,this%jne)+1 fld2(i,:) = fld(i,:)*float( icnt(this%ine,this%jne) )*float( i ) grd%tmp(this%ine,this%jne)=grd%tmp(this%ine,this%jne)+time_hash(this)*pos_hash(this)+log(this%mass) ichk4=ichk4+iberg ichk5=ichk5+iberg*i this=>this%next enddo ichk1=mpp_chksum( fld ) ichk2=mpp_chksum( fld2 ) ichk3=mpp_chksum( grd%tmp ) call mpp_sum( ichk4 ) call mpp_sum( ichk5 ) nbergs=count_bergs(bergs) if (nbergs.ne.sum(icnt(:,:))) then write(*,'("diamonds, bergs_chksum: ",2(a,i))') & '# bergs =', nbergs, ' sum(icnt) =',sum(icnt(:,:)) call error_mesg('diamonds, bergs_chksum:', 'mismatch in berg count!', FATAL) endif check_halo=.true. if (present(ignore_halo_violation)) then if (ignore_halo_violation) check_halo=.false. endif if (check_halo.and.nbergs.ne.sum(icnt(grd%isc:grd%iec, grd%jsc:grd%jec))) then write(*,'("diamonds, bergs_chksum: ",2(a,i))') & '# bergs =', nbergs, ' sum(icnt(comp_dom)) =',sum(icnt(:,:)) call error_mesg('diamonds, bergs_chksum:', 'mismatch in berg count on computational domain!', FATAL) endif call mpp_sum(nbergs) if (mpp_pe().eq.mpp_root_pe()) & write(*,'("diamonds, bergs_chksum: ",a18,6(x,a,"=",i))') & txt, 'chksum', ichk1, 'chksum2', ichk2, 'chksum3', ichk3, 'chksum4', ichk4, 'chksum5', ichk5, '#', nbergs deallocate( fld ) deallocate( fld2 ) deallocate( icnt ) if (debug) call count_out_of_order(bergs,txt) end subroutine bergs_chksum ! ############################################################################## integer function berg_chksum(berg ) ! Arguments type(iceberg), pointer :: berg ! Local variables real :: rtmp(28) integer :: itmp(28+3), i8=0, ichk1, ichk2, ichk3 integer :: i rtmp(:)=0. rtmp(1)=berg%lon rtmp(2)=berg%lat rtmp(3)=berg%uvel rtmp(4)=berg%vvel rtmp(5)=berg%mass rtmp(6)=berg%thickness rtmp(7)=berg%width rtmp(8)=berg%length rtmp(9)=berg%start_lon rtmp(10)=berg%start_lat rtmp(11)=berg%start_day rtmp(12)=berg%start_mass rtmp(13)=berg%mass_scaling rtmp(14)=berg%mass_of_bits rtmp(15)=berg%heat_density rtmp(16)=berg%xi rtmp(17)=berg%yj rtmp(19)=berg%uo rtmp(20)=berg%vo rtmp(21)=berg%ui rtmp(22)=berg%vi rtmp(23)=berg%ua rtmp(24)=berg%va rtmp(25)=berg%ssh_x rtmp(26)=berg%ssh_y rtmp(27)=berg%cn rtmp(28)=berg%hi itmp(1:28)=transfer(rtmp,i8) itmp(29)=berg%start_year itmp(30)=berg%ine itmp(31)=berg%jne ichk1=0; ichk2=0; ichk3=0 do i=1,28+3 ichk1=ichk1+itmp(i) ichk2=ichk2+itmp(i)*i ichk3=ichk3+itmp(i)*i*i enddo berg_chksum=ichk1+ichk2+ichk3 end function berg_chksum ! ############################################################################## end module