<p><b>dwj07@fsu.edu</b> 2012-10-15 14:33:08 -0600 (Mon, 15 Oct 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding a logical to control monthly forcing.<br>
Adding a lighter weight version of cullLoops.<br>
Changing double precision to real*8 because gfortran doesn't like having the types mismatched.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/basin_pbc_mrp/src/basin.F
===================================================================
--- branches/ocean_projects/basin_pbc_mrp/src/basin.F 2012-10-15 20:30:41 UTC (rev 2209)
+++ branches/ocean_projects/basin_pbc_mrp/src/basin.F 2012-10-15 20:33:08 UTC (rev 2210)
@@ -63,7 +63,7 @@
real(kind=4), allocatable, dimension(:,:,:) :: TAUX_MONTHLY, TAUY_MONTHLY
real, dimension(40) :: dz
-integer, parameter :: nMonths = 12
+integer :: nMonths = 1
! Step 1: Set the number of Vertical levels
integer, parameter :: nVertLevelsMOD = 40
@@ -91,6 +91,7 @@
logical, parameter :: eliminate_inland_seas=.true.
logical, parameter :: l_woce = .true.
logical, parameter :: write_OpenDX_flag = .true.
+logical, parameter :: monthly_forcing = .false.
! Set the number of top layers that are not allowed to have land, usually three.
integer, parameter :: top_layers_without_land = 3
@@ -148,10 +149,18 @@
integer :: iMonth
character(len=80) :: fileNameT, fileNameS, fileNameU
-fileNameT = 'TS/woce_t_ann.3600x2431x42interp.r4.nc'
-fileNameS = 'TS/woce_s_ann.3600x2431x42interp.r4.nc'
-fileNameU = 'TS/ws.old_ncep_1958-2000avg.interp3600x2431.nc'
+if(monthly_forcing) then
+ nMonths = 12
+else
+ nMonths = 1
+end if
+
+
+fileNameT = 'TS/annual/woce_t_ann.3600x2431x42interp.r4.nc'
+fileNameS = 'TS/annual/woce_s_ann.3600x2431x42interp.r4.nc'
+fileNameU = 'TS/annual/ws.old_ncep_1958-2000avg.interp3600x2431.nc'
+
! get to work
write(6,*) ' starting'
write(6,*)
@@ -205,118 +214,119 @@
call read_TS_fields(t_lon, t_lat, depth_t, TEMP, SALT, TAUX, TAUY)
call read_TS_finalize()
-allocate(SST_MONTHLY(nlon,nlat,12), SSS_MONTHLY(nlon,nlat,12))
-allocate(TAUX_MONTHLY(nlon,nlat,12), TAUY_MONTHLY(nlon,nlat,12))
-SST_MONTHLY=0; SSS_MONTHLY=0; TAUX_MONTHLY=0; TAUY_MONTHLY=0
+if(monthly_forcing) then
+ allocate(SST_MONTHLY(nlon,nlat,nMonths), SSS_MONTHLY(nlon,nlat,nMonths))
+ allocate(TAUX_MONTHLY(nlon,nlat,nMonths), TAUY_MONTHLY(nlon,nlat,nMonths))
+ SST_MONTHLY=0; SSS_MONTHLY=0; TAUX_MONTHLY=0; TAUY_MONTHLY=0
+ iMonth=1
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.01.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly01.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.01.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=2
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.02.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly02.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.02.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=3
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.03.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly03.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.03.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=4
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.04.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly04.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.04.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=5
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.05.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly05.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.05.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=6
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.06.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly06.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.06.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=7
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.07.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly07.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.07.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=8
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.08.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly08.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.08.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=9
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.09.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly09.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.09.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=10
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.10.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly10.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.10.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=11
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.11.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly11.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.11.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+
+ iMonth=12
+ fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.12.interp3600x2431.nc'
+ fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly12.interp3600x2431.nc'
+ fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.12.interp3600x2431.nc'
+ call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
+ write(6,*) nlon,nlat,ndepth
+ call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
+ call read_MONTHLY_finalize()
+end if
-iMonth=1
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.01.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly01.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.01.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=2
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.02.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly02.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.02.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=3
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.03.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly03.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.03.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=4
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.04.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly04.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.04.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=5
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.05.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly05.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.05.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=6
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.06.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly06.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.06.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=7
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.07.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly07.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.07.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=8
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.08.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly08.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.08.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=9
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.09.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly09.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.09.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=10
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.10.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly10.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.10.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=11
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.11.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly11.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.11.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
-iMonth=12
-fileNameT='TS/SST.shf.normal_year+Hurrell.monthly.12.interp3600x2431.nc'
-fileNameS='TS/SSS.sfwf.CORE_SSS+precip.monthly12.interp3600x2431.nc'
-fileNameU='TS/TAUIJ.ws.old_ncep_1958-2000avg.mon.12.interp3600x2431.nc'
-call read_MONTHLY_init(nlon, nlat, ndepth,fileNameT, fileNameS, fileNameU)
-write(6,*) nlon,nlat,ndepth
-call read_MONTHLY_fields(SST_MONTHLY(:,:,iMonth), SSS_MONTHLY(:,:,iMonth), TAUX_MONTHLY(:,:,iMonth), TAUY_MONTHLY(:,:,iMonth))
-call read_MONTHLY_finalize()
-
do k=1,ndepth
ndata = 0; temp_t=0; temp_s=0
do j=1,nlat
@@ -671,28 +681,32 @@
u_srcNew(1,iEdge) = ux*q(1) + uy*q(2) + uz*q(3)
endif
- do iMonth=1,12
- ulon = TAUX_MONTHLY(ix,iy,iMonth)
- ulat = TAUY_MONTHLY(ix,iy,iMonth)
-! if(abs(ulon).gt.1.0.or.abs(ulat).gt.1.0) then
-! ulon=0.0
-! ulat=0.0
-! endif
- call transform_from_lonlat_to_xyz(xin,yin,zin,ulon,ulat,ux,uy,uz)
- if(boundaryEdgeNew(1,iEdge).eq.1) then
- windStressMonthlyNew(iMonth,iEdge) = 0.0
+ if(monthly_forcing) then
+ do iMonth=1,nMonths
+ ulon = TAUX_MONTHLY(ix,iy,iMonth)
+ ulat = TAUY_MONTHLY(ix,iy,iMonth)
+ ! if(abs(ulon).gt.1.0.or.abs(ulat).gt.1.0) then
+ ! ulon=0.0
+ ! ulat=0.0
+ ! endif
+ call transform_from_lonlat_to_xyz(xin,yin,zin,ulon,ulat,ux,uy,uz)
+ if(boundaryEdgeNew(1,iEdge).eq.1) then
+ windStressMonthlyNew(iMonth,iEdge) = 0.0
+ else
+ iCell1 = cellsOnEdgeNew(1,iEdge)
+ iCell2 = cellsOnEdgeNew(2,iEdge)
+ p(1) = xCellNew(iCell1); p(2) = yCellNew(iCell1); p(3) = zCellNew(iCell1)
+ q(1) = xCellNew(iCell2); q(2) = yCellNew(iCell2); q(3) = zCellNew(iCell2)
+ q = q - p
+ call unit_vector_in_3space(q)
+ ! repeat
+ windStressMonthlyNew(iMonth,iEdge) = ux*q(1) + uy*q(2) + uz*q(3)
+ ! windStressMonthlyNew(iMonth,iEdge) = u_srcNew(1,iEdge)
+ endif
+ enddo
else
- iCell1 = cellsOnEdgeNew(1,iEdge)
- iCell2 = cellsOnEdgeNew(2,iEdge)
- p(1) = xCellNew(iCell1); p(2) = yCellNew(iCell1); p(3) = zCellNew(iCell1)
- q(1) = xCellNew(iCell2); q(2) = yCellNew(iCell2); q(3) = zCellNew(iCell2)
- q = q - p
- call unit_vector_in_3space(q)
- ! repeat
- windStressMonthlyNew(iMonth,iEdge) = ux*q(1) + uy*q(2) + uz*q(3)
- ! windStressMonthlyNew(iMonth,iEdge) = u_srcNew(1,iEdge)
- endif
- enddo
+ windStressMonthlyNew(:,:) = 0.0
+ end if
enddo
@@ -791,78 +805,83 @@
temperatureRestoreNew(:) = temperatureNew(1,1,:)
salinityRestoreNew(:) = salinityNew(1,1,:)
-do iMonth=1,12
-iNoData = 0
-do iCell=1,nCellsNew
- if(mod(iCell,100).eq.0) write(6,*) 'l_woce t and s RESTORE',iCell
- rlon = lonCellNew(iCell)/dtr
- rlat = latCellNew(iCell)/dtr
- ix = nint(rlon/0.1 - 0.05) + nlon + 1
- ix = mod(ix,nlon)+1
- iy = nlat
- do j=1,nlat
- if(t_lat(j).gt.rlat) then
- iy = j
- exit
- endif
- enddo ! j
- k=1
- ndata = 0; temp_t = 0; temp_s = 0
- do i=-15,15
- ixt = ix + 8*i
- if(ixt.lt.1) then
- ixt = ixt + nlon
- elseif(ixt.gt.nlon) then
- ixt = ixt - nlon
+if(monthly_forcing) then
+ do iMonth=1,nMonths
+ iNoData = 0
+ do iCell=1,nCellsNew
+ if(mod(iCell,100).eq.0) write(6,*) 'l_woce t and s RESTORE',iCell
+ rlon = lonCellNew(iCell)/dtr
+ rlat = latCellNew(iCell)/dtr
+ ix = nint(rlon/0.1 - 0.05) + nlon + 1
+ ix = mod(ix,nlon)+1
+ iy = nlat
+ do j=1,nlat
+ if(t_lat(j).gt.rlat) then
+ iy = j
+ exit
+ endif
+ enddo ! j
+ k=1
+ ndata = 0; temp_t = 0; temp_s = 0
+ do i=-15,15
+ ixt = ix + 8*i
+ if(ixt.lt.1) then
+ ixt = ixt + nlon
+ elseif(ixt.gt.nlon) then
+ ixt = ixt - nlon
+ endif
+ do j=-15,15
+ iyt = iy + 8*j
+ flag_lat = .true.
+ if(iyt.lt.1.or.iyt.gt.nlat) then
+ iyt = 1
+ flag_lat = .false.
+ endif
+ if(SST_MONTHLY(ixt,iyt,iMonth).gt.-10.0.and.flag_lat) then
+ ndata = ndata + 1
+ temp_t = temp_t + SST_MONTHLY(ixt,iyt,iMonth)
+ temp_s = temp_s + SSS_MONTHLY(ixt,iyt,iMonth)
+ endif
+ enddo !j
+ enddo !i
+
+ if(ndata.gt.0) then
+ temperatureRestoreMonthlyNew(iMonth,iCell) = temp_t / float(ndata)
+ salinityRestoreMonthlyNew(iMonth,iCell) = temp_s / float(ndata)
+ else
+ temperatureRestoreMonthlyNew(iMonth,iCell) = temperatureNew(1,1,iCell)
+ salinityRestoreMonthlyNew(iMonth,iCell) = salinityNew(1,1,iCell)
endif
- do j=-15,15
- iyt = iy + 8*j
- flag_lat = .true.
- if(iyt.lt.1.or.iyt.gt.nlat) then
- iyt = 1
- flag_lat = .false.
- endif
- if(SST_MONTHLY(ixt,iyt,iMonth).gt.-10.0.and.flag_lat) then
- ndata = ndata + 1
- temp_t = temp_t + SST_MONTHLY(ixt,iyt,iMonth)
- temp_s = temp_s + SSS_MONTHLY(ixt,iyt,iMonth)
- endif
- enddo !j
- enddo !i
-
- if(ndata.gt.0) then
- temperatureRestoreMonthlyNew(iMonth,iCell) = temp_t / float(ndata)
- salinityRestoreMonthlyNew(iMonth,iCell) = temp_s / float(ndata)
- else
- temperatureRestoreMonthlyNew(iMonth,iCell) = temperatureNew(1,1,iCell)
- salinityRestoreMonthlyNew(iMonth,iCell) = salinityNew(1,1,iCell)
- endif
-
+
+ enddo ! iCell
+ enddo ! iMonth
+
+ ! do a couple of smoothing passes
+ do iter=1,5
+ do iCell=1,nCellsNew
+ k=1
+ ndata=1
+ temp_t = temperatureRestoreMonthlyNew(iMonth,iCell)
+ temp_s = salinityRestoreMonthlyNew(iMonth,iCell)
+ do j=1,nEdgesOnCellNew(iCell)
+ jCell = cellsOnCellNew(j,iCell)
+ if(jCell.gt.0) then
+ if(maxLevelCellNew(jCell).ge.k) then
+ temp_t = temp_t + temperatureRestoreMonthlyNew(iMonth,iCell)
+ temp_s = temp_s + salinityRestoreMonthlyNew(iMonth,iCell)
+ ndata = ndata + 1
+ endif
+ endif
+ enddo ! j
+ temperatureRestoreMonthlyNew(iMonth,iCell) = temp_t / ndata
+ salinityRestoreMonthlyNew(iMonth,iCell) = temp_s / ndata
enddo ! iCell
-enddo ! iMonth
+ enddo ! iter
+else
+ temperatureRestoreMonthlyNew(:,:) = 0.0
+ salinityRestoreMonthlyNew(:,:) = 0.0
+end if
-! do a couple of smoothing passes
-do iter=1,5
-do iCell=1,nCellsNew
- k=1
- ndata=1
- temp_t = temperatureRestoreMonthlyNew(iMonth,iCell)
- temp_s = salinityRestoreMonthlyNew(iMonth,iCell)
- do j=1,nEdgesOnCellNew(iCell)
- jCell = cellsOnCellNew(j,iCell)
- if(jCell.gt.0) then
- if(maxLevelCellNew(jCell).ge.k) then
- temp_t = temp_t + temperatureRestoreMonthlyNew(iMonth,iCell)
- temp_s = temp_s + salinityRestoreMonthlyNew(iMonth,iCell)
- ndata = ndata + 1
- endif
- endif
- enddo ! j
- temperatureRestoreMonthlyNew(iMonth,iCell) = temp_t / ndata
- salinityRestoreMonthlyNew(iMonth,iCell) = temp_s / ndata
-enddo ! iCell
-enddo ! iter
-
endif ! l_woce
!repeat
Modified: branches/ocean_projects/basin_pbc_mrp/src/module_cullLoops.F
===================================================================
--- branches/ocean_projects/basin_pbc_mrp/src/module_cullLoops.F 2012-10-15 20:30:41 UTC (rev 2209)
+++ branches/ocean_projects/basin_pbc_mrp/src/module_cullLoops.F 2012-10-15 20:33:08 UTC (rev 2210)
@@ -1,13 +1,13 @@
module cullLoops
- public :: eliminateLoops
+ public :: eliminateLoops
- contains
+ contains
- subroutine eliminateLoops(nCells,nEdges,nVertices,maxEdges,vertexDegree, &
- nEdgesOnCell, cellsOnCell, verticesOnEdge, cellsOnVertex, edgesOnCell, lonCell, latCell, &
- xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex, &
- KMT)
+ subroutine eliminateLoops(nCells,nEdges,nVertices,maxEdges,vertexDegree, &
+ nEdgesOnCell, cellsOnCell, verticesOnEdge, cellsOnVertex, edgesOnCell, lonCell, latCell, &
+ xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex, &
+ KMT)
implicit none
@@ -30,274 +30,55 @@
logical :: connected, atBoundary, moveSouth, moveEast, atGrenwich
real :: lat, rlat, rlon, rCenter(3), s(3), t(3), q(3), rCross, mylon, mylat, pi
- pi = 4.0*atan(1.0)
+ integer, dimension(:), pointer :: cellStack
+ integer, dimension(:), pointer :: oceanMask
+ integer :: iCellStart, nStack, addedCells
+ real :: latStart, lonStart
- ! we loop over all cells and count the number of edges in the loop containing iCell
- ! there is no coupling between iCell, so this can be inside an openMP directive
- iCellMask(:) = 0
- moveSouth = .true.
- do iCell=1,nCells
+ write(6,*) 'Culling inland seas.....'
- write(6,*) 'working on : ',iCell, KMT(iCell)
-
- ! skip over land cells
- if(KMT(iCell).eq.0) then
- write(6,*) ' skipping : ', iCell
- cycle
- endif
+ allocate(cellStack(nCells/2))
+ allocate(oceanMask(nCells))
- ! the working cell will be jCell, so set jCell=iCell to start
- jCell = iCell
- ! write(6,*) 'setting jCell: ', jCell
+ oceanMask = 0
+ addedCells = 0
- atBoundary=.false. ! are we at a boundary?
- lCell = -1 ! when at a boundary, what is the index of the land cell or our right?
- oCell = -1 ! when at a boundary, what is the index of the ocean cell to our left?
+ iCellStart = maxloc(kmt, dim=1)
- do while (.not.atBoundary)
+ write(6,*) 'Starting index. ', iCellStart
+ write(6,*) 'lat, lon: ', latCell(iCellStart), lonCell(iCellStart)
+ write(6,*) 'Starting kmt: ', kmt(iCellStart)
- ! check to see if any edges of jCell are along the boundary
- do i=1,nEdgesOnCell(jCell)
- kCell = cellsOnCell(i,jCell)
- if(KMT(kCell).eq.0) then
- lCell = kCell ! this is a land cell
- oCell = jCell ! this is an ocean cell
- atBoundary = .true.
- write(6,*) ' found boundary : ',lCell,oCell
- endif
- enddo
+ nStack = 1
+ cellStack(nStack) = iCellStart
+ oceanMask(iCellStart) = 1
+ addedCells = 1
- ! choose the next cell to be the one with the most southern latitude
- ! this jCell will only be used if atBoundary=.false., thus jCell must be an ocean cell
- if(moveSouth) then
- rlat = 10000.0
- mylat = latCell(jCell)
- do i=1,nEdgesOnCell(jCell)
- kCell = cellsOnCell(i,jCell)
- if(latCell(kCell).lt.rlat) then
- rlat = latCell(kCell)
- iSave = kCell
- endif
- enddo
- jCell = iSave
- endif
+ do while(nStack > 0)
+ oCell = cellStack(nStack)
+ nStack = nStack - 1
+ write(6,*) ' Working on cell ', oCell, addedCells, nStack
- ! if(moveSouth.and..not.atBoundary) write(6,*) ' pushing on to the south : ', jCell
-
- enddo ! .not.atBoundary
+ do i = 1, nEdgesOnCell(oCell)
+ iCell = cellsOnCell(i, oCell)
- ! OK, we hit a boundary ..... trace out the full loop in the CCW direction
- write(6,*) ' OK we hit a boundary, let us trace out this loop '
- write(6,*) ' ocean cell ', oCell, KMT(oCell)
- write(6,*) ' land cell ', lCell, KMT(lCell)
+ if(kmt(iCell) > 0 .and. oceanMask(iCell) == 0) then
+ nStack = nStack + 1
+ cellStack(nStack) = iCell
+ oceanMask(iCell) = 1
+ addedCells = addedCells + 1
+ end if
+ end do
+ end do
- ! start the counter at 1
- iEdgeCounter = 1
- edgeList(:) = 0
+ where(oceanMask == 0) kmt(:) = 0
- ! find the shared edge where we are starting and save the starting edge index
- iSharedEdge = sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
- iStartEdge = iSharedEdge
- edgeList(iEdgeCounter) = iSharedEdge
+ write(6,*) addedCells, ' total cells have been in the stack.'
+ write(6,*) 'Done culling inland seas.....'
- connected = .false.
- LeftTurns = 0; RightTurns = 0
- do while (.not.connected)
+ deallocate(cellStack)
+ deallocate(oceanMask)
- call moveAhead(xCell,yCell,zCell,xVertex,yVertex,zVertex, &
- oCell,lCell,iSharedEdge,nCells,nEdges,nVertices,maxEdges,vertexDegree, &
- verticesOnEdge,cellsOnVertex,iCellAhead)
+ end subroutine eliminateLoops
- ! if the cell ahead is ocean, then the boundary is shared between lCell and iCellAhead
- ! if the cell ahead is land, then the boundary is shared between oCell and iCellAhead
- if(KMT(iCellAhead).gt.0) then
- oCell = iCellAhead
- RightTurns = RightTurns + 1
- ! write(6,*) ' the cell ahead is ocean, will turn right ', lCell, oCell
- else
- lCell = iCellAhead
- LeftTurns = LeftTurns + 1
- ! write(6,*) ' the cell ahead is land, will turn left ', lCell, oCell
- endif
- iSharedEdge = sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
-
- ! check to see if we are where we started
- if(iSharedEdge.eq.iStartEdge) then
- connected=.true.
- write(6,*) ' we are back to the where we started '
- else
- iEdgeCounter=iEdgeCounter+1
- edgeList(iEdgeCounter) = iSharedEdge
- endif
-
- enddo ! .not.connected
-
- ! OK, we now have a loop .... but did we circle an inland see or a land mass?
- rCenter(:) = 0.0
- do iEdge=1,iEdgeCounter
- rCenter(1) = rCenter(1) + xEdge(edgeList(iEdge))/iEdgeCounter
- rCenter(2) = rCenter(2) + yEdge(edgeList(iEdge))/iEdgeCounter
- rCenter(3) = rCenter(3) + zEdge(edgeList(iEdge))/iEdgeCounter
- enddo
- rCenter(:) = rCenter(:) / sqrt ( rCenter(1)**2 + rCenter(2)**2 + rCenter(3)**2 )
-
- rCross = 0.0
- do iEdge=1,iEdgeCounter-1
- t(1) = xEdge(edgeList(iEdge+1)) - xEdge(edgeList(iEdge))
- t(2) = yEdge(edgeList(iEdge+1)) - yEdge(edgeList(iEdge))
- t(3) = zEdge(edgeList(iEdge+1)) - zEdge(edgeList(iEdge))
- s(1) = rCenter(1) - xEdge(edgeList(iEdge))
- s(2) = rCenter(2) - yEdge(edgeList(iEdge))
- s(3) = rCenter(3) - zEdge(edgeList(iEdge))
- t(:) = t(:) / sqrt( t(1)**2 + t(2)**2 + t(3)**2 )
- s(:) = s(:) / sqrt( s(1)**2 + s(2)**2 + s(3)**2 )
- call cross_product_in_R3(t,s,q)
- rCross = rCross + q(1)*rCenter(1) + q(2)*rCenter(2) + q(3)*rCenter(3)
- enddo
-
- write(6,*)
- write(6,*) ' edges and cross ', iEdgeCounter, rCross, LeftTurns, RightTurns
- write(6,*)
-
- if(LeftTurns-RightTurns.gt.0.and.rCross.gt.0.0) then
- iCellMask(iCell) = 1
- write(50,11) iCell, lonCell(iCell), latCell(iCell)
- 11 format(i8,2f12.4)
- endif
-
- enddo ! iCell
-
- ! cull all inland seas
- do iSweep=1,nCells/10
- write(6,*) iSweep
- do iCell=1,nCells
- if(iCellMask(iCell).eq.1) then
- do i=1,nEdgesOnCell(iCell)
- kCell=cellsOnCell(i,iCell)
- if(KMT(kCell).gt.0) iCellMask(kCell)=1
- enddo
- endif
- enddo
- enddo
-
- write(6,*) ' total cells culled ', sum(iCellMask)
- where(iCellMask(:).eq.1) KMT(:)=0
-
- end subroutine eliminateLoops
-
-
- subroutine moveAhead(xCell,yCell,zCell,xVertex,yVertex,zVertex, &
- oCell,lCell,iSharedEdge,nCells,nEdges,nVertices,maxEdges,vertexDegree, &
- verticesOnEdge,cellsOnVertex,iCellAhead)
- implicit none
- integer, intent(in) :: oCell, lCell,iSharedEdge, nCells, nEdges, nVertices, maxEdges, vertexDegree
- integer, intent(in) :: verticesOnEdge(2,nEdges), cellsOnVertex(vertexDegree,nVertices)
- real, intent(in), dimension(nCells) :: xCell, yCell, zCell
- real, intent(in), dimension(nVertices) :: xVertex, yVertex, zVertex
- integer, intent(out) :: iCellAhead
- integer :: iVertex1,iVertex2, i, kCell
- real :: v1(3), v2(3), v3(3), ocean(3), land(3), d1, d2, cross1(3), cross2(3)
-
- ! solution assumes a CCW ordering of cellsOnVertex
- ! the cell ahead will be connected to the vertex that lists cellsOnVertex with lCell following oCell
-
- ! the vertex moving in the CCW direction has to be one of the two vertices connected to iSharedEdge
- iVertex1 = verticesOnEdge(1,iSharedEdge)
- iVertex2 = verticesOnEdge(2,iSharedEdge)
- !write(6,*) ' iVertex1, iVertex2 ', iVertex1, iVertex2
- !write(6,*) cellsOnVertex(:,iVertex1)
- !write(6,*) cellsOnVertex(:,iVertex2)
-
- v1(1)=xVertex(iVertex1)
- v1(2)=yVertex(iVertex1)
- v1(3)=zVertex(iVertex1)
-
- v2(1)=xVertex(iVertex2)
- v2(2)=yVertex(iVertex2)
- v2(3)=zVertex(iVertex2)
-
- ocean(1) = xCell(oCell)
- ocean(2) = yCell(oCell)
- ocean(3) = zCell(oCell)
-
- land(1) = xCell(lCell)
- land(2) = yCell(lCell)
- land(3) = zCell(lCell)
-
- v1 = v1 - ocean
- v2 = v2 - ocean
- v3 = land - ocean
-
- ocean(:) = ocean(:) / sqrt( ocean(1)**2 + ocean(2)**2 + ocean(3)**2)
- land(:) = land(:) / sqrt( land(1)**2 + land(2)**2 + land(3)**2)
- v1(:) = v1(:) / sqrt( v1(1)**2 + v1(2)**2 + v1(3)**2)
- v2(:) = v2(:) / sqrt( v2(1)**2 + v2(2)**2 + v2(3)**2)
- v3(:) = v3(:) / sqrt( v3(1)**2 + v3(2)**2 + v3(3)**2)
-
- call cross_product_in_R3(v3,v1,cross1)
- call cross_product_in_R3(v3,v2,cross2)
-
- d1 = ocean(1)*cross1(1) + ocean(2)*cross1(2) + ocean(3)*cross1(3)
- d2 = ocean(1)*cross2(1) + ocean(2)*cross2(2) + ocean(3)*cross2(3)
-
- kCell = 0
-
- ! write(6,*) lCell, oCell
- ! write(6,*) cellsOnVertex(:,iVertex1), d1
- ! write(6,*) cellsOnVertex(:,iVertex2), d2
- ! write(6,*) v1(:)
- ! write(6,*) v2(:)
-
- if(d1.gt.d2) then
- do i=1,vertexDegree
- kCell=cellsOnVertex(i,iVertex1)
- if(kCell.ne.oCell.and.kCell.ne.lCell) iCellAhead=kCell
- enddo
- endif
-
- if(d2.gt.d1) then
- do i=1,vertexDegree
- kCell=cellsOnVertex(i,iVertex2)
- if(kCell.ne.oCell.and.kCell.ne.lCell) iCellAhead=kCell
- enddo
- endif
-
- end subroutine moveAhead
-
-
-
- function sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
- implicit none
- integer, intent(in) :: oCell, lCell, nCells, maxEdges, nEdgesOnCell(nCells), edgesOnCell(maxEdges,nCells)
- integer :: sharedEdge
- integer :: i,j,iEdge,jEdge
-
- sharedEdge=-1
- do i=1,nEdgesOnCell(oCell)
- iEdge = edgesOnCell(i,oCell)
- do j=1,nEdgesOnCell(lCell)
- jEdge = edgesOnCell(j,lCell)
- if(iEdge.eq.jEdge) then
- sharedEdge = jEdge
- exit
- endif
- enddo
- enddo
-
- if(SharedEdge.eq.-1) then
- write(6,*) ' problem with SharedEdge ',oCell,lCell
- stop
- endif
-
- end function sharedEdge
-
- subroutine cross_product_in_R3(p_1,p_2,p_out)
- real , intent(in) :: p_1 (3), p_2 (3)
- real , intent(out) :: p_out (3)
- p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
- p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
- p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
- end subroutine cross_product_in_R3
-
-
end module cullLoops
Modified: branches/ocean_projects/basin_pbc_mrp/src/module_write_netcdf.F
===================================================================
--- branches/ocean_projects/basin_pbc_mrp/src/module_write_netcdf.F 2012-10-15 20:30:41 UTC (rev 2209)
+++ branches/ocean_projects/basin_pbc_mrp/src/module_write_netcdf.F 2012-10-15 20:33:08 UTC (rev 2210)
@@ -103,7 +103,7 @@
integer, intent(in) :: vertexDegree
integer, intent(in) :: nMonths
character (len=16) :: on_a_sphere
- double precision :: sphere_radius
+ real*8 :: sphere_radius
integer :: nferr
</font>
</pre>