particleTracking/readOcean.f90 at master · mlicer/particleTracking · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
! --------------------------------- read ocean T S RHO etc ------------

subroutine readOCEAN
	use ptmModule
	use netcdf
	use utilities
	implicit none
	include 'netcdf.inc'

	character*255 :: filename
	real, allocatable :: tmp(:,:,:,:)

	!      integer  :: nlon, nlat, nLevels
	integer  :: file_id,t_id,s_id, rho_id, hdiff_id, vdiff_id     ! variable dimensions
	integer  :: i,j,t, hour
	real :: missValnum
	integer  :: status, hour_t, hour_s,hour_rho, hour_hdiff, hour_vdiff
	character*10 ::missVal


	allocate(tmp(ncgrid%nLon,ncgrid%nLat,ncgrid%nLevels,setup%modelRuntimeInHours))

	! 3D version:
	!    allocate(sea%T(ncgrid%nLon,ncgrid%nLat,ncgrid%nLevels,setup%numHoursOfTracking))
	!    allocate(sea%S(ncgrid%nLon,ncgrid%nLat,ncgrid%nLevels,setup%numHoursOfTracking))
	!    allocate(sea%RHO(ncgrid%nLon,ncgrid%nLat,ncgrid%nLevels,setup%numHoursOfTracking))
	!    allocate(sea%HorizontalDiffusivity(ncgrid%nLon,ncgrid%nLat,ncgrid%nLevels,setup%numHoursOfTracking))
	!    allocate(sea%VerticalDiffusivity(ncgrid%nLon,ncgrid%nLat,ncgrid%nLevels,setup%numHoursOfTracking))

	allocate(sea%T(ncgrid%nLon,ncgrid%nLat,1,setup%numHoursOfTracking))
	allocate(sea%S(ncgrid%nLon,ncgrid%nLat,1,setup%numHoursOfTracking))
	allocate(sea%RHO(ncgrid%nLon,ncgrid%nLat,1,setup%numHoursOfTracking))
	allocate(sea%HorizontalDiffusivity(ncgrid%nLon,ncgrid%nLat,1,setup%numHoursOfTracking))
	allocate(sea%VerticalDiffusivity(ncgrid%nLon,ncgrid%nLat,1,setup%numHoursOfTracking))

	! set correct indexes for reading variables
	hour_t = 0
	hour_hdiff = 0
	hour_vdiff = 0
	!    hour_s = 0
	!    hour_rho = 0

	do i = 1,setup%numOfExistingFiles

	 filename=setup%existingFiles(i)

	 print *,"READING OCEAN STATE FROM: ",trim(adjustl(filename))

	 call handle_err(NF90_OPEN(trim(adjustl(filename)), NF90_NOWRITE, file_id))
	 call handle_err(NF90_INQ_VARID(file_id, 'temperature', t_id))
	 ! reading horizontal diffusivity ID
	 ! Note that in our POM horizontal kinematic viscosity = horizontal diffusivity since TPRNI = 1.
	 call handle_err(NF90_INQ_VARID(file_id, 'horizontal_kinematic_viscosity', hdiff_id))
	 ! reading vertical diffusivity ID
	 call handle_err(NF90_INQ_VARID(file_id, 'vertical_diffusivity', vdiff_id))
	!    call handle_err(NF90_INQ_VARID(file_id, 'salinity', s_id))

	! reading temperature and diffusivities:
	 call handle_err(NF90_GET_VAR(file_id, t_id, tmp))
	 call handle_err(NF90_GET_ATT(file_id, t_id, "missing_value", missVal))


	 call rstr2num(missVal,missValnum)

	if (.not. setup%performBacktracking) then
		do t = setup%hourReadStart(i),setup%hourReadEnd(i)
		    hour_t=hour_t+1
		    !sea%T(:,:,:,hour_t)=tmp(:,:,:,t)
		    sea%T(:,:,1,hour_t)=tmp(:,:,setup%level,t)
		enddo
	else
		do t = setup%hourReadEnd(i),setup%hourReadStart(i),-1
		    hour_t=hour_t+1
		    print *, "HOUR OF DAY, ARRAY INDEX:", t, hour_t
	!            sea%T(:,:,:,hour_t)=tmp(:,:,:,t)
		    sea%T(:,:,1,hour_t)=tmp(:,:,setup%level,t)
		enddo
	endif

	! reading horizontal diffusivity
	call handle_err(NF90_GET_VAR(file_id, hdiff_id, tmp))

	if (.not. setup%performBacktracking) then
		do t = setup%hourReadStart(i),setup%hourReadEnd(i)
		    hour_hdiff=hour_hdiff+1
		    !sea%HorizontalDiffusivity(:,:,:,hour_hdiff)=tmp(:,:,:,t)
		    sea%HorizontalDiffusivity(:,:,1,hour_hdiff)=tmp(:,:,setup%level,t)
		enddo
	else
		do t = setup%hourReadEnd(i),setup%hourReadStart(i),-1
		    hour_hdiff=hour_hdiff+1
	!            sea%HorizontalDiffusivity(:,:,:,hour_hdiff)=tmp(:,:,:,t)
		    sea%HorizontalDiffusivity(:,:,1,hour_hdiff)=tmp(:,:,setup%level,t)
		enddo
	endif

	! reading vertical diffusivity
	call handle_err(NF90_GET_VAR(file_id, vdiff_id, tmp))
	if (.not. setup%performBacktracking) then
		do t = setup%hourReadStart(i),setup%hourReadEnd(i)
		    hour_vdiff=hour_vdiff+1
		    !sea%VerticalDiffusivity(:,:,:,hour_hdiff)=tmp(:,:,:,t)
		    sea%VerticalDiffusivity(:,:,1,hour_hdiff)=tmp(:,:,setup%level,t)
		enddo
	else
		do t = setup%hourReadEnd(i),setup%hourReadStart(i),-1
		    hour_vdiff=hour_vdiff+1
	!            sea%VerticalDiffusivity(:,:,:,hour_vdiff)=tmp(:,:,:,t)
		    sea%VerticalDiffusivity(:,:,1,hour_vdiff)=tmp(:,:,setup%level,t)
		enddo
	endif

	 ! close file
	 call handle_err(NF90_CLOSE(file_id))

	enddo


end subroutine readOCEAN