forked from nansencenter/enkf-topaz
-
Notifications
You must be signed in to change notification settings - Fork 0
/
m_get_mod_fld.F90
199 lines (172 loc) · 7.12 KB
/
m_get_mod_fld.F90
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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
module m_get_mod_fld
! KAL -- This routine reads one of the fields from the model, specified
! KAL -- by name, vertical level and time level
! KAL -- This routine is really only effective for the new restart files.
contains
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
subroutine get_mod_fld(fld,j,cfld,vlevel,tlevel,nx,ny)
#if defined (QMPI)
use qmpi
#else
use qmpi_fake
#endif
implicit none
integer, intent(in) :: nx,ny ! Grid dimension
integer, intent(in) :: j ! Ensemble member to read
real, dimension(nx,ny), intent(out) :: fld ! output fld
character(len=*), intent(in) :: cfld ! name of fld
integer, intent(in) :: vlevel ! vertical level
integer, intent(in) :: tlevel ! time level
integer reclICE
real*8, dimension(nx,ny) :: ficem,hicem,hsnwm,ticem,tsrfm
logical ex
character(len=*),parameter :: icefile='forecastICE.uf'
! KAL -- shortcut -- the analysis is for observation icec -- this little "if"
! means the analysis will only work for ice. Add a check though
if ((trim(cfld)/='icec' .and. trim(cfld)/='hice') .or. vlevel/=0 .or. tlevel/=1)then
if (master) print *,'get_mod_fld only works for icec for now '//trim(cfld)
stop
end if
!###################################################################
!####################### READ ICE MODEL #########################
!###################################################################
#if defined (ICE)
#warning "COMPILING WITH ICE"
inquire(exist=ex,file=icefile)
if (.not.ex) then
if (master) then
print *,icefile//' does not exist!'
print *,'(get_mod_fld)'
end if
stop
end if
inquire(iolength=reclICE)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
open(10,file=icefile,form='unformatted',access='direct',recl=reclICE,action='read')
read(10,rec=j)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
if (trim(cfld)=='icec') fld = ficem
if (trim(cfld)=='hice') fld = hicem
close(10)
#else
#warning "COMPILING WITHOUT ICE"
#endif
return
end subroutine get_mod_fld
! KAL - This is for the new file type
subroutine get_mod_fld_new(memfile,fld,iens,cfld0,vlevel,tlevel,nx,ny,Indfield)
use mod_raw_io
use m_get_mod_fld_nc
#if defined (QMPI)
use qmpi, only : qmpi_proc_num, master,stop_mpi
#else
use qmpi_fake
#endif
implicit none
integer, intent(in) :: nx,ny ! Grid dimension
integer, intent(in) :: iens ! Ensemble member to read
real, dimension(nx,ny), intent(out) :: fld ! output fld
character(len=*), intent(in) :: memfile! base name of input files
character(len=*), intent(in) :: cfld0 ! name of fld
integer, intent(in) :: vlevel ! vertical level
integer, intent(in) :: tlevel ! time level
integer, intent(in) :: Indfield ! index in the analysisfields
real*8, dimension(nx,ny) :: readfldr8
real*4, dimension(nx,ny) :: readfldr4
real*4:: amin, amax,spval
real :: bmin, bmax
integer :: indx
!-----------------------------------
character*80 :: cfld,icefile
cfld=trim(cfld0)
#if defined (HYCOM_CICE)
if (trim(cfld0) == 'SSH') then
cfld='srfhgt'
end if
#endif
#if defined (SINGLE_RESTART)
if (trim(cfld0) == 'icec') then
cfld='ficem'
elseif (trim(cfld0) == 'hice') then
cfld='hicem'
endif
#endif
!-----------------------------------
! Dette fordi is-variablane forelobig er paa gammalt format.
if (trim(cfld) /= 'icec' .and. trim(cfld) /= 'hice' .and. &
trim(cfld) /= 'hsnwm' .and. trim(cfld) /= 'aicen' .and. &
trim(cfld) /= 'vicen' .and. trim(cfld) /= 'vsnon' .and. &
trim(cfld) /= 'ticem' .and. trim(cfld) /= 'tsrfm' .and. &
trim(cfld) /= 'sicem' .and. trim(cfld) /= 'tsnom' .and. &
trim(cfld) /= 'ficem' .and. trim(cfld) /= 'hicem') then
! KAL - 1) f kva index som skal lesast finn vi fraa .b fil (header)
! working well for restart file
indx=0;
if (Indfield>=0) then
call rst_index_from_header(trim(memfile)//'.b', & ! filnavn utan extension
cfld , & ! felt som skal lesast fex saln,temp
vlevel, & ! vertikalnivaa
tlevel, & ! time level - kan vere 1 eller 2 - vi bruker 1 foreloepig
indx, & ! indexen som maa lesas fra data fila
bmin,bmax, & ! min/max - kan sjekkast mot det som er i datafila
.true. )
else
call daily_index_from_header(trim(memfile)//'.b', & ! filnavn utan extension
cfld , & ! felt som skal lesast fex saln,temp
vlevel, & ! vertikalnivaa
indx, & ! indexen som maa lesas fra data fila
bmin,bmax,iens) ! min/max - kan sjekkast mot det som er i datafila
endif
if (indx < 0) then
if (master) then
print *, 'ERROR: get_mod_fld_new() to read ', trim(memfile), '.b: "',&
trim(cfld), '" not found for ',indx
end if
stop
end if
! KAL -- les datafelt vi fann fraa header fila (indx)
spval=2**100
call READRAW(readfldr4 ,& ! Midlertidig felt som skal lesast
amin, amax ,& ! max/min fraa data (.a) fila
nx,ny ,& ! dimensjonar
.true.,spval ,& ! dette brukast for sette "no value" verdiar
trim(memfile)//'.a',& ! fil som skal lesast fraa
indx) ! index funne over
! Sjekk p at vi har lest rett - samanlign max/min fr filene
if (abs(amin-bmin).gt.abs(bmin)*1.e-4 .or. &
abs(bmax-amax).gt.abs(bmax)*1.e-4 ) then
print *,'Inconsistency between .a and .b files (m_get_mod_fld)'
print *,'.a : ',amin,amax
print *,'.b : ',bmin,bmax
print *,trim(cfld),vlevel,tlevel, indx
print *,'node ',qmpi_proc_num
call exit(1)
end if
fld=readfldr4
else ! fld = fice, hice
! Gammal rutine ja
#if defined(HYCOM_CICE)
indx=-1
icefile=trim(memfile)//'.nc'
call get_mod_fld_nc(trim(icefile), readfldr8, &
cfld, vlevel, indx, nx, ny)
if (indx==-1) then
icefile='ice_'//trim(memfile)//'.nc'
call get_mod_fld_nc(trim(icefile), readfldr8, &
cfld, vlevel, indx, nx, ny)
end if
if (indx < 0) then
if (master) then
print *, 'ERROR: get_mod_fld_new(): ', trim(icefile),&
trim(cfld), '" not found for ', vlevel
stop
end if
stop
! else
! print *, trim(icefile), ' '//trim(cfld), vlevel,Indfield
end if
#else
call get_mod_fld(readfldr8,iens,cfld,0,1,nx,ny)
#endif
fld=readfldr8
end if
end subroutine
end module m_get_mod_fld