Commit 9d839238 authored by ph290's avatar ph290
Browse files

removing an issue where it still seemed to have a spinup day run and included...

removing an issue where it still seemed to have a spinup day run and included in the outout at he start of each year
parent 3f3c6f2e
......@@ -1095,9 +1095,16 @@ newdepth = alldepth
!if(alldepth.lt.10) newN = 5
!if(alldepth.lt.5) newN = 3
!kb=0.0025
kb=0.0025
!All vertical levels set to 2m
if(alldepth.gt.0.0) then
newN=nint(alldepth/2.0)
newdepth=2.0*newN
else
newN=1
newdepth=0.0
endif
!set up equal depth intervals of 2m for sections:
if(imode.eq.2) then
......@@ -1113,7 +1120,7 @@ endif
!open(1,file='s12_m2_s2_n2_h.dat',status='old')
!open(1,file='s12_m2_s2_n2_h_sec.dat',status='old')
!open(1,file='s12_m2_s2_n2_h_tim.dat',status='old')
vismax=0.1; Nz_bg=1.0e-5; Kz_bg=1.0e-5; first_temp=5.0; lambda=0.1; heat_shade=0.012; bed_din=woa_nutrient; din_rate=10.0
vismax=0.1; Nz_bg=1.0e-5; Kz_bg=1.0e-5; first_temp=25.0; lambda=0.1; heat_shade=0.012; bed_din=woa_nutrient; din_rate=10.0
par_atten=0.1; par_percent=0.45
uamp(1)=0.4; uamp(2:5)=0.0; uphase(1:5)=0.0; vphase(1:5)=0.0
vamp(1:5)=0.0
......@@ -1326,10 +1333,10 @@ Subroutine grazing() ! Get phytoplankton grazing parameters
metfile_error=0
end if
do i=1,ndays
day_angle=(360.0*dble(i)/dble(ndays)) * (dble(ndays)/365.0)
! day_angle=(360.0*dble(i)/dble(ndays)) * (dble(ndays)/365.0)
!The change in the line above alows for simulations longer than 1 year, but assumes all years are 365 days long.
declination(i)=0.39637-22.9133*MyCOS(day_angle)+4.02543*MySIN(day_angle)-0.3872*&
MyCOS(2.0*day_angle)+0.052*MySIN(2.0*day_angle)
! declination(i)=0.39637-22.9133*MyCOS(day_angle)+4.02543*MySIN(day_angle)-0.3872*&
! MyCOS(2.0*day_angle)+0.052*MySIN(2.0*day_angle)
end do
end if
......@@ -1341,7 +1348,7 @@ Subroutine grazing() ! Get phytoplankton grazing parameters
cloud(ndays+1)=cloud(ndays)
rad_input(ndays+1)=rad_input(ndays)
lw_rad_down(ndays+1)= lw_rad_down(ndays)
declination(ndays+1)=declination(ndays)
! declination(ndays+1)=declination(ndays)
1303 return
end subroutine get_met
......@@ -1457,7 +1464,7 @@ N=newN; depth=newdepth; dz=newdz
! avoid shallow and deep ocean
if(depth.lt.5.0.or.depth.gt.150.0) then
iday = 0
iday = 1
strat_jul = -999.99
timeloop_init: do itime=0,ndays
write(6,fmt="(i4,2f9.3,4f8.2)")iday,lon,lat,depth,strat_jul,strat_jul,strat_jul
......@@ -1476,8 +1483,8 @@ N=newN; depth=newdepth; dz=newdz
!
! Set physical arrays to initial values.
!
velx_old(1:200)=0.0; vely_old(1:200)=0.0; temp_old(1:200)=first_temp; velx_new(1:200)=0.0
vely_new(1:200)=0.0; temp_new(1:200)=first_temp; Nz(1:200)=vismax/2.0; Kz(1:200)=vismax/2.0
velx_old(1:200)=0.0; vely_old(1:200)=0.0; temp_old(1:200)=first_temp; velx_new(1:200)=0.1
vely_new(1:200)=0.1; temp_new(1:200)=first_temp; Nz(1:200)=vismax/2.0; Kz(1:200)=vismax/2.0
tke(1:N) = 3.0e-6
eps(1:N) = 5.0e-10
rad_amp=radiation(1); radsum=0.0; wx=1.0; wy=1.0
......@@ -1598,7 +1605,7 @@ end do
!now removed because model starts from restart files
!
! starttime=-72.0*3600.0; idmet=1; iday=-3; julian_day=-3; imonth=1; montest=1
starttime=0.0; idmet=1; iday=0; julian_day=0; imonth=1; montest=1
starttime=0.0; idmet=1; iday=1; julian_day=0; imonth=1; montest=1
!
rad_sum=0.0
......@@ -1641,21 +1648,22 @@ timeloop: do itime=1,itotal ! <A NAME="START OF TIME LOOP">
!
! Calculate u-cubed
!
if(iday.gt.-2)then
! if(iday.gt.-2)then
! if(iday.gt.2)then
speed_mag(1:N)=dsqrt(velx_new(1:N)**2.0+vely_new(1:N)**2.0)
mean_speed=sum(speed_mag(1:N))/real(N)
u_cubed=u_cubed+mean_speed**3.
n_cubed=n_cubed+1
end if
! end if
!
! Update radiation
!
if(iday.gt.1)idmet=iday
! if(iday.gt.1)idmet=iday
idmet=iday
noon_time=dble(idmet-1)+0.5
day_time=(time/(24.0*3600.0))-noon_time
day_angle=(2.0*3.1415927)*day_time
sin_solar_elev=dcos(day_angle)*MyCOS(lat)*MyCOS(declination(idmet))+MySIN(lat)*MySIN(declination(idmet))
! day_angle=(2.0*3.1415927)*day_time
! sin_solar_elev=dcos(day_angle)*MyCOS(lat)*MyCOS(declination(idmet))+MySIN(lat)*MySIN(declination(idmet))
!now directly reading downwelling radiation in
rad0 = rad_input(idmet)
lw_rad_down0 = lw_rad_down(idmet)
......@@ -2115,108 +2123,108 @@ month_stress=(month_stress+sqrt(stressx**2.0+stressy**2.0))/real(montest)
!print*, ' lon lat depth log(h/u3) strat accumC '
! map output
if(imode.eq.1) then
!write(6,'(2f8.3,5f8.2)') lon,lat,depth,dlog10(depth/u3_mean),strat_jul,rad_sum/acount,month_net1+accumulated,bottom_phyto_biomass
! write(6,fmt="(i4,2f8.3,5f8.2)")iday,lon,lat,depth,temp_new(N),temp_new(1),x_new(N),x_new(1)/chl_carbon
!note 2f8.3,8f8.2 means that there are two columns with 8 spaces for digits,
!three numbers below the decimal place, then 8 columns with 8 spaces for
!and two numbers behind the decimal place
write(6,fmt="(i4,2f8.3)",advance="no")iday,lon,lat
if(include_depth_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")depth
end if
if(include_temp_surface_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")temp_new(N)
end if
if(include_temp_bottom_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")temp_new(1)
end if
if(include_chlorophyll_surface_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")x_new(N)
end if
if(include_phyto_biomass_surface_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")x_new(N)/chl_carbon
end if
if(include_phyto_biomass_bottom_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")x_new(1)/chl_carbon
end if
if(include_PAR_surface_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")rad_mean(N)
end if
if(include_PAR_bottom_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")rad_mean(1)
end if
if(include_windspeed_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")windspeed
end if
if(include_stressx_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")stressx
end if
if(include_stressy_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")stressy
end if
if(include_Etide_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")Etide
end if
if(include_Ewind_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")Ewind
end if
if(include_u_mean_surface_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")u_mean(N)
end if
if(include_u_mean_bottom_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")u_mean(1)
end if
if(include_grow1_mean_surface_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")grow1_mean(N)
end if
if(include_grow1_mean_bottom_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")grow1_mean(1)
end if
if(include_uptake1_mean_surface_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")uptake1_mean(N)
end if
if(include_uptake1_mean_bottom_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")uptake1_mean(1)
end if
if(include_tpn1_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")tpn1
end if
if(include_tpg1_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")tpg1
end if
if(include_speed3_output.eq.1) then
write(6,fmt="(1f8.2)",advance="no")speed3
end if
write(6,fmt="()")
endif
! ! map output
! if(imode.eq.1) then
! !write(6,'(2f8.3,5f8.2)') lon,lat,depth,dlog10(depth/u3_mean),strat_jul,rad_sum/acount,month_net1+accumulated,bottom_phyto_biomass
! ! write(6,fmt="(i4,2f8.3,5f8.2)")iday,lon,lat,depth,temp_new(N),temp_new(1),x_new(N),x_new(1)/chl_carbon
! !note 2f8.3,8f8.2 means that there are two columns with 8 spaces for digits,
! !three numbers below the decimal place, then 8 columns with 8 spaces for
! !and two numbers behind the decimal place
!
! write(6,fmt="(i4,2f8.3)",advance="no")iday,lon,lat
!
!
! if(include_depth_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")depth
! end if
!
! if(include_temp_surface_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")temp_new(N)
! end if
!
! if(include_temp_bottom_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")temp_new(1)
! end if
!
! if(include_chlorophyll_surface_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")x_new(N)
! end if
!
! if(include_phyto_biomass_surface_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")x_new(N)/chl_carbon
! end if
!
! if(include_phyto_biomass_bottom_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")x_new(1)/chl_carbon
! end if
!
! if(include_PAR_surface_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")rad_mean(N)
! end if
!
! if(include_PAR_bottom_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")rad_mean(1)
! end if
!
! if(include_windspeed_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")windspeed
! end if
!
! if(include_stressx_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")stressx
! end if
!
! if(include_stressy_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")stressy
! end if
!
! if(include_Etide_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")Etide
! end if
!
! if(include_Ewind_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")Ewind
! end if
!
! if(include_u_mean_surface_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")u_mean(N)
! end if
!
! if(include_u_mean_bottom_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")u_mean(1)
! end if
!
! if(include_grow1_mean_surface_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")grow1_mean(N)
! end if
!
! if(include_grow1_mean_bottom_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")grow1_mean(1)
! end if
!
! if(include_uptake1_mean_surface_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")uptake1_mean(N)
! end if
!
! if(include_uptake1_mean_bottom_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")uptake1_mean(1)
! end if
!
! if(include_tpn1_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")tpn1
! end if
!
! if(include_tpg1_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")tpg1
! end if
!
! if(include_speed3_output.eq.1) then
! write(6,fmt="(1f8.2)",advance="no")speed3
! end if
!
! write(6,fmt="()")
!
! endif
! time series output
if(imode.eq.3) then
......@@ -2435,7 +2443,10 @@ spec_hum1=0.62*sat_vap/(airP(idmet)-0.38*sat_vap)
spec_hum2=0.62*vap/(airP(idmet)-0.38*vap)
!hl=0.985*5.67d-8*(surf_temp**4.0)*(0.39-0.05*vap**0.5)*(1.0-0.6d-4*cloud(idmet)**2.0) !original codes Longwave
!Net longwave radiation flux=IR_surface_emissivity_1_for_black_body(Stefan_Boltzmann_constant * T**4.8 - downwelling_longwave) Note T should really be skin temp
hl=0.985*5.67d-8*(surf_temp**4.0)-lw_rad_down0 !Longwave
!hl=0.985*5.67d-8*(surf_temp**4.0)-(lw_rad_down0) !Longwave
hl=0.99*5.67d-8*(surf_temp**4.0)-(lw_rad_down0*0.75) !Longwave The 0.95 is
!playing around to account for lw being absorbed and loast in teh skin rather
!than top meters
sk=1.45d-3*1.3*1004.0*wind_speed(idmet)*(temp_old(N)-airT(idmet)) !Sensible
ske=1.5d-3*1.3*wind_speed(idmet)*(spec_hum1-spec_hum2)*(2.5d6-2.3d3*temp_old(N)) !Evaporation
rad_out=-(hl+sk+ske)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment