Added A68 test back into repository and driver. Tested it with new nm…

…l option 'use_broken_bonds_for_substep_contact=.true.', which yielded about a 15% speedup compared to previously

driver/icebergs_driver.F90

@@ -62,11 +62,15 @@ program icebergs_driver
 logical :: grounding_test=.false.
 logical :: big_grounding_test=.false.
 logical :: dem_test_contact=.false.
+logical :: a68_test=.false.
 logical :: vert_int_ocean_vel=.false.
 logical :: tau_is_velocity=.true.
+logical :: reverse_a68_forcings=.false.
 logical :: sin_test=.false.
 logical :: fl_test=.false.
 character(len=1000) :: data_dir = 'data/'
+integer :: transient_a68_data_start_ind=0 !<set to .gt. 0 to use the hourly JRA-55 wind data for the A68
+!! test. The values indicate the time index to start from in the hourly data
 integer :: halo=1 !< Width of halo in parallel decomposition
 real :: ibdt=3600.0
 real :: ibuo=0.0
@@ -87,8 +91,8 @@ program icebergs_driver
 integer :: write_time_inc=1
 namelist /icebergs_driver_nml/ debug, ni, nj, halo, ibhrs, ibdt, ibuo, ibvo, nmax, &
   saverestart,ibui,ibvi,collision_test,chaotic_test,grounding_test,&
-  gridres,write_time_inc,bump_depth, sst, data_dir, vert_int_ocean_vel,&
-  sin_test,fl_test,tau_is_velocity,ibua,ibuy,&
+  gridres,write_time_inc,bump_depth, sst, a68_test, data_dir, vert_int_ocean_vel,&
+  reverse_a68_forcings,sin_test,fl_test,tau_is_velocity,transient_a68_data_start_ind,ibua,ibuy,&
   Old_wrong_Rearth_and_SSH,no_wind,REarth,big_grounding_test,dem_test_contact !mom_Rearth
 ! For loops
 integer :: isc !< Start of i-index for computational domain (used for loops)
@@ -123,6 +127,8 @@ program icebergs_driver
 real, dimension(:,:), allocatable :: uo,vo !<zonal and meridonal ocean velocities (m/s)
 real, dimension(:,:), allocatable :: ui,vi !<zonal and meridonal ice velocities (m/s)
 real, dimension(:,:), allocatable :: tauxa,tauya !< Zonal and meridonal wind stress (Pa)
+real, dimension(:,:,:), allocatable :: tauxa_hr,tauya_hr !< Hourly jra-55 zonal and meridonal wind velocity (m/s) for A68 test
+real, dimension(:,:,:), allocatable :: uo_hr,vo_hr !<Hourly OSCAR zonal and meridonal ocean surface velocities (m/s) for A68 test
 real, dimension(:,:), allocatable :: ssh !< Effective sea-surface height (m)
 real, dimension(:,:,:), allocatable :: ssh_hr !< Effective sea-surface height (m, hourly)
 real, dimension(:,:), allocatable :: sstemp !< Sea-surface temperature (C or K)
@@ -147,6 +153,12 @@ program icebergs_driver
 read(nmunit, nml=icebergs_driver_nml)
 call close_file(nmunit)
 
+if (a68_test) then
+  call a68_dims(data_dir,ni,nj)
+  if (mpp_pe()==0) print *,'ni,nj',ni,nj
+  iflags=0
+endif
+
 if (min(ni,nj)<1) call error_mesg('icebergs_driver:', &
   'Both ni,nj must be positive in icebergs_driver_nml namelist', FATAL)
 
@@ -241,91 +253,140 @@ program icebergs_driver
 ! ustar_berg = !Friction velocity on base of bergs (m/s)
 ! area_berg = !Area of bergs (m2)
 
-do j = jsd, jed
-  do i = isd, ied
-    lon(i,j) = gridres*real(i)
-    lat(i,j) = gridres*real(j)
-    dx(i,j) = gridres
-    dy(i,j) = gridres
-    area(i,j) = gridres*gridres
-    wet(i,j) = 1.0
-    cos_rot(i,j) = 1.0
-    sin_rot(i,j) = 0.0
-    depth(i,j) = 1000.
-  enddo
-enddo
+if (a68_test) then
+  if (Old_wrong_Rearth_and_SSH) REarth=6360000.
+  call mpp_sync()
+  call a68_prep(data_dir,mpp_domain,lon,lat,dx,dy,area,depth,uo,vo,tauxa,tauya,ssh,vert_int_ocean_vel,tau_is_velocity,Old_wrong_Rearth_and_SSH,REarth)
+  !mom_Rearth)
+  call mpp_sync()
+  if (no_wind) then
+    tauxa=0.0; tauya=0.0
+  endif
+  !where (lon<(-37.6+360.))
+  !  depth=depth+1000.
+  !end where
+  wet=1.0
+  cos_rot=1.0
+  sin_rot=0.0
+  if (reverse_a68_forcings) then
+    uo=-uo; vo=-vo
+    tauxa=-tauxa; tauya=-tauya
+  endif
 
-if (chaotic_test) then
-  !assign land cells at the N and S portions of the domain.
-  !input.nml: set coastal drift parameter to prevent grounding
-  do j=jsd,jed; do i=isd,ied
-    if (lat(i,j)<=0.0+2*gridres .or. lat(i,j)>=1000.e3-2*gridres) wet(i,j)=0.0
-  enddo;enddo
-endif
+  if (transient_a68_data_start_ind>0) then
+    if (.not. (ibdt==3600.0 .or. ibdt==1800.0)) then
+      call error_mesg('icebergs_driver:', &
+      'in order to use transient_a68_data_start_ind>0, timestep must be 30 min or 1 hr', FATAL)
+    endif
+    if (.not. tau_is_velocity) call error_mesg('icebergs_driver:', &
+      'tau_is_velocity must be true if transient_a68_data_start_ind>-1', FATAL)
+    allocate( tauxa_hr(isd:ied,jsd:jed,739) )
+    allocate( tauya_hr(isd:ied,jsd:jed,739) )
+    allocate( uo_hr(isd:ied,jsd:jed,744) )
+    allocate( vo_hr(isd:ied,jsd:jed,744) )
+    allocate( ssh_hr(isd:ied,jsd:jed,744) )
+    call a68_prep_3d(data_dir,mpp_domain,tauxa_hr,tauya_hr,uo_hr,vo_hr,ssh_hr)
+    if (reverse_a68_forcings) then
+      tauxa_hr=-tauxa_hr; tauya_hr=-tauya_hr
+      uo_hr=-uo_hr; vo_hr=-vo_hr
+      ssh_hr=-ssh_hr
+      dit=-1
+    else
+      dit=1
+    endif
+    if (no_wind) then
+      tauxa_hr=0.0; tauya_hr=0.0
+    endif
+  endif
+else
+  do j = jsd, jed
+    do i = isd, ied
+      lon(i,j) = gridres*real(i)
+      lat(i,j) = gridres*real(j)
+      dx(i,j) = gridres
+      dy(i,j) = gridres
+      area(i,j) = gridres*gridres
+      wet(i,j) = 1.0
+      cos_rot(i,j) = 1.0
+      sin_rot(i,j) = 0.0
+      depth(i,j) = 1000.
+    enddo
+  enddo
 
-if (big_grounding_test) then
-  lat(:,:)=lat(:,:)-0.45 !-25000.1
-  lon(:,:)=lon(:,:)-0.45
-  !assign land cells at the N and S portions of the domain.
-  !input.nml: set coastal drift parameter to prevent grounding
-  do j=jsd,jed; do i=isd,ied
-    if (lat(i,j)<=(-5.e3) .or. lat(i,j)>=(220.e3)) wet(i,j)=0.0
-  enddo;enddo
+  if (chaotic_test) then
+    !assign land cells at the N and S portions of the domain.
+    !input.nml: set coastal drift parameter to prevent grounding
+    do j=jsd,jed; do i=isd,ied
+      if (lat(i,j)<=0.0+2*gridres .or. lat(i,j)>=1000.e3-2*gridres) wet(i,j)=0.0
+    enddo;enddo
+  endif
 
-  !introduce gaussian bump to bathymetry
-  a = 1000.0-bump_depth !max height(m)
-  c = 5e3; !controls width
-  bx = 63.e3; by = 60.e3 !50.e3
-  do j=jsd,jed; do i=isd,ied
-    xc=lon(i,j)-(gridres/2.); yc=lat(i,j)-(gridres/2.) !define at cell centers
-    depth(i,j)=a*exp(-((xc-bx)*(xc-bx)/(2.*c*c)+(yc-by)*(yc-by)/(2.*c*c)))
-  enddo; enddo
-  print *,'min, max depth',minval(depth(:,:)),maxval(depth(:,:))
-  depth = 1000.0-depth
-endif
+  if (big_grounding_test) then
+    lat(:,:)=lat(:,:)-0.45 !-25000.1
+    lon(:,:)=lon(:,:)-0.45
+    !assign land cells at the N and S portions of the domain.
+    !input.nml: set coastal drift parameter to prevent grounding
+    do j=jsd,jed; do i=isd,ied
+      if (lat(i,j)<=(-5.e3) .or. lat(i,j)>=(220.e3)) wet(i,j)=0.0
+    enddo;enddo
+
+    !introduce gaussian bump to bathymetry
+    a = 1000.0-bump_depth !max height(m)
+    c = 5e3; !controls width
+    bx = 63.e3; by = 60.e3 !50.e3
+    do j=jsd,jed; do i=isd,ied
+      xc=lon(i,j)-(gridres/2.); yc=lat(i,j)-(gridres/2.) !define at cell centers
+      depth(i,j)=a*exp(-((xc-bx)*(xc-bx)/(2.*c*c)+(yc-by)*(yc-by)/(2.*c*c)))
+    enddo; enddo
+    print *,'min, max depth',minval(depth(:,:)),maxval(depth(:,:))
+    depth = 1000.0-depth
+  endif
 
-if (dem_test_contact) then
-  lat(:,:)=lat(:,:)-0.45 !-25000.1
-  lon(:,:)=lon(:,:)-0.45
-  !assign land cells at the N and S portions of the domain.
-  !input.nml: set coastal drift parameter to prevent grounding
-  do j=jsd,jed; do i=isd,ied
-    if (lat(i,j)<=(-5.e3) .or. lat(i,j)>=(220.e3)) wet(i,j)=0.0
-  enddo;enddo
+   if (dem_test_contact) then
+    lat(:,:)=lat(:,:)-0.45 !-25000.1
+    lon(:,:)=lon(:,:)-0.45
+    !assign land cells at the N and S portions of the domain.
+    !input.nml: set coastal drift parameter to prevent grounding
+    do j=jsd,jed; do i=isd,ied
+      if (lat(i,j)<=(-5.e3) .or. lat(i,j)>=(220.e3)) wet(i,j)=0.0
+    enddo;enddo
 
-  mid=50.e3
+    mid=50.e3
 
-  do j=jsd,jed; do i=isd,ied
-    if (lon(i,j)>65 .or. lat(i,j)==mid) then
-      vo(i,j)=0.0
-    else
-      if (lat(i,j)>mid) then
-        vo(i,j)=-ibvo
+    do j=jsd,jed; do i=isd,ied
+      if (lon(i,j)>65 .or. lat(i,j)==mid) then
+        vo(i,j)=0.0
       else
-        vo(i,j)=ibvo
+        if (lat(i,j)>mid) then
+          vo(i,j)=-ibvo
+        else
+          vo(i,j)=ibvo
+        endif
       endif
-    endif
   enddo;enddo
-endif
 
-if (grounding_test) then
-  lat(:,:)=lat(:,:)-15000.0
-  !assign land cells at the N and S portions of the domain.
-  !input.nml: set coastal drift parameter to prevent grounding
-  do j=jsd,jed; do i=isd,ied
-    if (lat(i,j)<=0.0+2*gridres .or. lat(i,j)>=30.e3-2*gridres) wet(i,j)=0.0
-  enddo;enddo
 
-  !introduce gaussian bump to bathymetry
-  a = 1000.0-bump_depth !max height(m)
-  c = 2.5e3; !controls width
-  bx = 15.e3; by = 15.e3
-  do j=jsd,jed; do i=isd,ied
-    xc=lon(i,j)-(gridres/2.); yc=lat(i,j)-(gridres/2.) !define at cell centers
-    depth(i,j)=a*exp(-((xc-bx)*(xc-bx)/(2.*c*c)+(yc-by)*(yc-by)/(2.*c*c)))
-  enddo; enddo
-  print *,'min, max depth',minval(depth(:,:)),maxval(depth(:,:))
-  depth = 1000.0-depth
+  endif
+
+  if (grounding_test) then
+    lat(:,:)=lat(:,:)-15000.0
+    !assign land cells at the N and S portions of the domain.
+    !input.nml: set coastal drift parameter to prevent grounding
+    do j=jsd,jed; do i=isd,ied
+      if (lat(i,j)<=0.0+2*gridres .or. lat(i,j)>=30.e3-2*gridres) wet(i,j)=0.0
+    enddo;enddo
+
+    !introduce gaussian bump to bathymetry
+    a = 1000.0-bump_depth !max height(m)
+    c = 2.5e3; !controls width
+    bx = 15.e3; by = 15.e3
+    do j=jsd,jed; do i=isd,ied
+      xc=lon(i,j)-(gridres/2.); yc=lat(i,j)-(gridres/2.) !define at cell centers
+      depth(i,j)=a*exp(-((xc-bx)*(xc-bx)/(2.*c*c)+(yc-by)*(yc-by)/(2.*c*c)))
+    enddo; enddo
+    print *,'min, max depth',minval(depth(:,:)),maxval(depth(:,:))
+    depth = 1000.0-depth
+  endif
 endif
 
 if (sin_test) vo = vo * sin(lat*2*pi/(20.0*gridres) + gridres)
@@ -336,7 +397,11 @@ program icebergs_driver
 
 call set_calendar_type(THIRTY_DAY_MONTHS)
 
-time = set_date(1,1,1,0,0,0)
+if (transient_a68_data_start_ind==355) then
+  time = set_date(1,1,7,5,0,0)
+else
+  time = set_date(1,1,1,0,0,0)
+endif
 dt = ibdt
 
 call icebergs_init(bergs, &
@@ -395,6 +460,7 @@ program icebergs_driver
 endif
 
 ns = 1 !timestep counter
+if (transient_a68_data_start_ind>0) ns2 = 1
 
 ! Time_end = increment_date(Time, years, months, days, hours, minutes, seconds)
 Time_end = increment_date(Time,0,0,0,ibhrs,0,0)
@@ -415,6 +481,23 @@ program icebergs_driver
     write(*,*) ''
   end if
 
+  if (transient_a68_data_start_ind>0) then
+
+    if ((ibdt==3600.0) .or. (mod(ns2,1.)==0)) then
+      tauxa=tauxa_hr(:,:,transient_a68_data_start_ind+dit*(int(ns2)-1))
+      tauya=tauya_hr(:,:,transient_a68_data_start_ind+dit*(int(ns2)-1))
+      uo=uo_hr(:,:,transient_a68_data_start_ind+dit*(int(ns2)-1))
+      vo=vo_hr(:,:,transient_a68_data_start_ind+dit*(int(ns2)-1))
+      ssh=ssh_hr(:,:,transient_a68_data_start_ind+dit*(int(ns2)-1))
+    else
+      tauxa=0.5*(tauxa+tauxa_hr(:,:,transient_a68_data_start_ind+dit*(int(ceiling(ns2))-1)))
+      tauya=0.5*(tauya+tauya_hr(:,:,transient_a68_data_start_ind+dit*(int(ceiling(ns2))-1)))
+      uo=0.5*(uo+uo_hr(:,:,transient_a68_data_start_ind+dit*(int(ceiling(ns2))-1)))
+      vo=0.5*(vo+vo_hr(:,:,transient_a68_data_start_ind+dit*(int(ceiling(ns2))-1)))
+      ssh=ssh_hr(:,:,transient_a68_data_start_ind+dit*(int(ns2)-1))
+    endif
+  endif
+
   ! The main driver the steps updates icebergs
   call icebergs_run(bergs, time, calving(isc:iec,jsc:jec), &
                     uo(isc-1:iec+1,jsc-1:jec+1), vo(isc-1:iec+1,jsc-1:jec+1), &
@@ -430,8 +513,41 @@ program icebergs_driver
 
   Time = Time + real_to_time_type(dt)
   ns = ns+1
+
+  if (transient_a68_data_start_ind>0) then
+    if (ibdt==3600.0) then
+      ns2=ns2+1
+    elseif (ibdt==1800.0) then
+      ns2=ns2+0.5
+    endif
+  endif
 enddo
 
+
+if (a68_test) then
+  call cpu_time(time_finish)
+  call get_date(Time, iyr, imon, iday, ihr, imin, isec)
+  write(*,*) ''
+  write(*,'(a)') '-------------------------------------------'
+  write(*,'(a,i5)')' Saving after timestep   ',ns-1
+  print *,' ns2   ',ns2
+  if (mod(ns2,1.)==0) then
+    print *,'transient_a68_data_start_ind+dit*(int(ns2)-1)',&
+      transient_a68_data_start_ind+dit*(int(ns2)-1)
+  else
+    print *,'transient_a68_data_start_ind+dit*(int(ceiling(ns2))-1))',&
+      transient_a68_data_start_ind+dit*(int(ceiling(ns2))-1)
+  endif
+  write(*,'(a)')' Restart time is:'
+  write(*,'(a,i5,a,i5,a,i5)')' year',iyr,'  month ',imon,'  day   ',iday
+  write(*,'(a,i5,a,i5,a,i5)')' hour',ihr,'  minute',imin,'  second',isec
+  write(*,*) ''
+  write(*,*) 'clock-time elapsed: ', (time_finish - time_begin)/60., 'minutes'
+  write(*,*) 'clock-time per day: ', (time_finish - time_begin)/(ns-1)*(60*60*24)/ibdt
+  write(*,'(a)') '-------------------------------------------'
+  write(*,*) ''
+endif
+
 if (saverestart) call icebergs_save_restart(bergs)
 
 call cpu_time(time_finish)