Skip to content

Commit 66f9d86

Browse files
jjokellajjokella
committed
Merge branch 'master' into dev-snow-da-up-to-date-old
2 parents 47c6e93 + bae512a commit 66f9d86

File tree

12 files changed

+426
-394
lines changed

12 files changed

+426
-394
lines changed

‎bldsva/intf_DA/pdaf/framework/Makefile‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -145,7 +145,7 @@ $(PROG) : $(LIBMODEL) libpdaf-d.a \
145145
$(MODULES) $(MOD_ASSIM) $(OBJ_MODEL_PDAF) $(OBJ_PDAF_INT) $(OBJ_PDAF_USER)
146146
$(PREP_C) $(LD) $(OPT_LNK) -o $@ \
147147
$(MODULES) $(MOD_ASSIM) $(OBJ_MODEL_PDAF) $(OBJ_PDAF_INT) $(OBJ_PDAF_USER) \
148-
-L$(BASEDIR)/lib -lpdaf-d $(LINK_LIBS) $(LIBMODEL) $(LIBS)
148+
-L$(BASEDIR)/lib -lpdaf-d $(LIBMODEL) $(LIBS) $(LINK_LIBS)
149149

150150
######################################################
151151

‎bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -120,7 +120,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
120120
use ColumnType, only : col
121121
! use GetGlobalValuesMod, only: GetGlobalWrite
122122
! use clm_varcon, only: nameg
123-
use enkf_clm_mod, only: col_index_hydr_act
123+
use enkf_clm_mod, only: state_clm2pdaf_p
124124
use enkf_clm_mod, only: clmstatevec_only_active
125125
use enkf_clm_mod, only: clmstatevec_max_layer
126126
#else
@@ -788,7 +788,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
788788
endif
789789
enddo
790790

791-
if(obs_interp_switch) then
791+
if(obs_interp_switch.eq.1) then
792792
! loop over all obs and save the indices of the nearest grid
793793
! points to array obs_interp_indices_p and save the distance
794794
! weights to array obs_interp_weights_p (later normalized)
@@ -939,7 +939,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
939939
! vector assembling.
940940
if(clmstatevec_allcol.eq.1) then
941941
#ifdef CLMFIVE
942-
if(clmstatevec_only_active) then
942+
if(clmstatevec_only_active.eq.1) then
943943

944944
! Error if observation deeper than clmstatevec_max_layer
945945
if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then
@@ -951,7 +951,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
951951
print *, "clmstatevec_max_layer=", clmstatevec_max_layer
952952
call abort_parallel()
953953
end if
954-
obs_index_p(cnt) = col_index_hydr_act(c,clmobs_layer(i))
954+
obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i))
955955
else
956956
#endif
957957
obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1))
@@ -978,7 +978,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
978978
end do
979979
end do
980980

981-
if(obs_interp_switch) then
981+
if(obs_interp_switch.eq.1) then
982982
! loop over all obs and save the indices of the nearest grid
983983
! points to array obs_interp_indices_p and save the distance
984984
! weights to array obs_interp_weights_p (later normalized)

‎bldsva/intf_DA/pdaf/framework/init_dim_obs_l_pdaf.F90‎

Lines changed: 16 additions & 114 deletions
Original file line numberDiff line numberDiff line change
@@ -129,75 +129,14 @@ SUBROUTINE init_dim_obs_l_pdaf(domain_p, step, dim_obs_f, dim_obs_l)
129129
#endif
130130

131131
! Count observations within cradius
132-
#ifdef CLMSA
133-
obsind = 0
134-
obsdist = 0.0
135-
dim_obs_l = 0
136-
if(point_obs.eq.0) then
137-
max_var_id = MAXVAL(var_id_obs_nc(:,:))
138-
allocate(log_var_id(max_var_id))
139-
log_var_id(:) = .TRUE.
140132

141-
do m = 1, dim_nx
142-
do k = 1, dim_ny
143-
i = (m-1)* dim_ny + k
144-
do j = 1, max_var_id
145-
if(log_var_id(j) .and. var_id_obs_nc(k,m) == j) then
146-
dx = abs(longxy_obs(i) - longxy(domain_p))
147-
dy = abs(latixy_obs(i) - latixy(domain_p))
148-
dist = sqrt(real(dx)**2 + real(dy)**2)
149-
!obsdist(i) = dist
150-
if(dist == 0) then
151-
dim_obs_l = dim_obs_l + 1
152-
obsind(i) = 1
153-
log_var_id(j) = .FALSE.
154-
obsdist(i) = dist
155-
EXIT
156-
end if
157-
end if
158-
end do
159-
enddo
160-
enddo
161-
162-
do m = 1, dim_nx
163-
do k = 1, dim_ny
164-
i = (m-1)* dim_ny + k
165-
do j = 1, max_var_id
166-
if(log_var_id(j) .and. var_id_obs_nc(k,m) == j) then
167-
dx = abs(longxy_obs(i) - longxy(domain_p))
168-
dy = abs(latixy_obs(i) - latixy(domain_p))
169-
dist = sqrt(real(dx)**2 + real(dy)**2)
170-
!obsdist(i) = dist
171-
if (dist <= real(cradius)) then
172-
dim_obs_l = dim_obs_l + 1
173-
obsind(i) = 1
174-
log_var_id(j) = .FALSE.
175-
dx = abs(lon_var_id(j) - longxy(domain_p))
176-
dy = abs(lat_var_id(j) - latixy(domain_p))
177-
obsdist(i) = sqrt(real(dx)**2 + real(dy)**2)
178-
end if
179-
end if
180-
end do
181-
enddo
182-
enddo
183-
else
184-
do i = 1,dim_obs
185-
dx = abs(longxy_obs(i) - longxy(domain_p))
186-
dy = abs(latixy_obs(i) - latixy(domain_p))
187-
dist = sqrt(real(dx)**2 + real(dy)**2)
188-
obsdist(i) = dist
189-
if (dist <= real(cradius)) then
190-
dim_obs_l = dim_obs_l + 1
191-
obsind(i) = 1
192-
end if
193-
end do
194-
end if
195-
#else
196-
#ifdef PARFLOW_STAND_ALONE
133+
#ifndef CLMSA
134+
#ifndef OBS_ONLY_CLM
197135
obsind = 0
198136
obsdist = 0.0
199137
dim_obs_l = 0
200138
if(point_obs.eq.0) then
139+
if(model == tag_model_parflow) THEN
201140
max_var_id = MAXVAL(var_id_obs_nc(:,:))
202141
allocate(log_var_id(max_var_id))
203142
log_var_id(:) = .TRUE.
@@ -228,7 +167,9 @@ SUBROUTINE init_dim_obs_l_pdaf(domain_p, step, dim_obs_f, dim_obs_l)
228167
end do
229168
end do
230169
end do
170+
end if
231171
else
172+
if(model == tag_model_parflow) THEN
232173
do i = 1,dim_obs
233174
dx = abs(x_idx_obs_nc(i) - int(xcoord_fortran(domain_p_coord))-1)
234175
dy = abs(y_idx_obs_nc(i) - int(ycoord_fortran(domain_p_coord))-1)
@@ -239,45 +180,21 @@ SUBROUTINE init_dim_obs_l_pdaf(domain_p, step, dim_obs_f, dim_obs_l)
239180
obsind(i) = 1
240181
end if
241182
end do
183+
endif
242184
endif
243-
#else
185+
#endif
186+
#endif
187+
188+
#ifndef PARFLOW_STAND_ALONE
189+
#ifndef OBS_ONLY_PARFLOW
244190
obsind = 0
245191
obsdist = 0.0
246192
dim_obs_l = 0
247193
if(point_obs.eq.0) then
248-
if(model == tag_model_parflow) THEN
249194
max_var_id = MAXVAL(var_id_obs_nc(:,:))
250195
allocate(log_var_id(max_var_id))
251196
log_var_id(:) = .TRUE.
252197

253-
do m = 1, dim_nx
254-
do k = 1, dim_ny
255-
i = (m-1)* dim_ny + k
256-
do j = 1, max_var_id
257-
if(log_var_id(j) .and. var_id_obs_nc(k,m) == j) then
258-
dx = abs(x_idx_obs_nc(i) - int(xcoord_fortran(domain_p_coord))-1)
259-
dy = abs(y_idx_obs_nc(i) - int(ycoord_fortran(domain_p_coord))-1)
260-
dist = sqrt(real(dx)**2 + real(dy)**2)
261-
!obsdist(i) = dist
262-
if (dist <= real(cradius) .AND. dist > 0) then
263-
dim_obs_l = dim_obs_l + 1
264-
obsind(i) = 1
265-
log_var_id(j) = .FALSE.
266-
dx = abs(ix_var_id(j) - int(xcoord_fortran(domain_p_coord))-1)
267-
dy = abs(iy_var_id(j) - int(ycoord_fortran(domain_p_coord))-1)
268-
obsdist(i) = sqrt(real(dx)**2 + real(dy)**2)
269-
else if(dist == 0) then
270-
dim_obs_l = dim_obs_l + 1
271-
obsind(i) = 1
272-
log_var_id(j) = .FALSE.
273-
obsdist(i) = dist
274-
end if
275-
end if
276-
end do
277-
end do
278-
end do
279-
end if
280-
281198
if(model == tag_model_clm) THEN
282199
do m = 1, dim_nx
283200
do k = 1, dim_ny
@@ -322,22 +239,7 @@ SUBROUTINE init_dim_obs_l_pdaf(domain_p, step, dim_obs_f, dim_obs_l)
322239
enddo
323240
enddo
324241
end if
325-
! print *,'mype_filter dim_obs_l ', mype_filter, dim_obs_l
326-
else
327-
if(model == tag_model_parflow) THEN
328-
do i = 1,dim_obs
329-
dx = abs(x_idx_obs_nc(i) - int(xcoord_fortran(domain_p_coord))-1)
330-
dy = abs(y_idx_obs_nc(i) - int(ycoord_fortran(domain_p_coord))-1)
331-
dist = sqrt(real(dx)**2 + real(dy)**2)
332-
obsdist(i) = dist
333-
if (dist <= real(cradius)) then
334-
dim_obs_l = dim_obs_l + 1
335-
obsind(i) = 1
336-
end if
337-
end do
338-
endif
339-
340-
#ifndef OBS_ONLY_PARFLOW
242+
else
341243
if(model == tag_model_clm) THEN
342244
do i = 1,dim_obs
343245
dx = abs(longxy_obs(i) - longxy(domain_p))
@@ -348,13 +250,13 @@ SUBROUTINE init_dim_obs_l_pdaf(domain_p, step, dim_obs_f, dim_obs_l)
348250
dim_obs_l = dim_obs_l + 1
349251
obsind(i) = 1
350252
end if
351-
end do
253+
end do
352254
end if
353-
#endif
354-
endif
255+
end if
355256
#endif
356257
#endif
357-
258+
!------------------------------------------------------------------------
259+
358260
! kuw: allocate and determine local observation index and distance
359261
!#ifndef CLMSA
360262
! Initialize index array for local observations in full observed vector

‎bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -116,7 +116,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
116116
use ColumnType, only : col
117117
! use GetGlobalValuesMod, only: GetGlobalWrite
118118
! use clm_varcon, only: nameg
119-
use enkf_clm_mod, only: col_index_hydr_act
119+
use enkf_clm_mod, only: state_clm2pdaf_p
120120
use enkf_clm_mod, only: clmstatevec_only_active
121121
use enkf_clm_mod, only: clmstatevec_max_layer
122122
#else
@@ -781,7 +781,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
781781
endif
782782
enddo
783783

784-
if(obs_interp_switch) then
784+
if(obs_interp_switch.eq.1) then
785785
! loop over all obs and save the indices of the nearest grid
786786
! points to array obs_interp_indices_p and save the distance
787787
! weights to array obs_interp_weights_p (later normalized)
@@ -932,7 +932,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
932932
! vector assembling.
933933
if(clmstatevec_allcol.eq.1) then
934934
#ifdef CLMFIVE
935-
if(clmstatevec_only_active) then
935+
if(clmstatevec_only_active.eq.1) then
936936

937937
! Error if observation deeper than clmstatevec_max_layer
938938
if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then
@@ -944,7 +944,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
944944
print *, "clmstatevec_max_layer=", clmstatevec_max_layer
945945
call abort_parallel()
946946
end if
947-
obs_index_p(cnt) = col_index_hydr_act(c,clmobs_layer(i))
947+
obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i))
948948
else
949949
#endif
950950
obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1))
@@ -971,7 +971,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
971971
end do
972972
end do
973973

974-
if(obs_interp_switch) then
974+
if(obs_interp_switch.eq.1) then
975975
! loop over all obs and save the indices of the nearest grid
976976
! points to array obs_interp_indices_p and save the distance
977977
! weights to array obs_interp_weights_p (later normalized)

‎bldsva/intf_DA/pdaf/framework/localize_covar_pdaf.F90‎

Lines changed: 44 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -98,6 +98,7 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH)
9898
real(r8), pointer :: lon(:)
9999
real(r8), pointer :: lat(:)
100100
integer, pointer :: mycgridcell(:) !Pointer for CLM3.5/CLM5.0 col->gridcell index arrays
101+
REAL :: yhalf
101102
#endif
102103
INTEGER :: icoord
103104

@@ -229,7 +230,23 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH)
229230
dx = abs(clmobs_lon(obs_pdaf2nc(j)) - lon(mycgridcell(state_pdaf2clm_c_p(i))))
230231
dy = abs(clmobs_lat(obs_pdaf2nc(j)) - lat(mycgridcell(state_pdaf2clm_c_p(i))))
231232

232-
distance = sqrt(real(dx)**2 + real(dy)**2)
233+
! Check for longitude differences that may yield differences
234+
! larger than 180deg depending on conventions. Example:
235+
! crossing the prime meridian (lon=0deg), when convention is
236+
! all-positive lons
237+
IF (dx > 180.0) THEN
238+
dx = 360.0 - dx
239+
END IF
240+
241+
! Intermediate latitude
242+
yhalf = ( clmobs_lat(obs_pdaf2nc(j)) + lat(mycgridcell(state_pdaf2clm_c_p(i))) ) / 2.0
243+
244+
! Latitude-dependent factor for longitude difference
245+
dx = dx * cos(yhalf * 3.14159265358979323846 / 180.0)
246+
247+
! Factor ca. 111km comes from R*pi/180, where R is earth
248+
! radius and pi/180 is because we input lat/lon in degrees
249+
distance = 111.19492664455873 * sqrt(real(dx)**2 + real(dy)**2)
233250

234251
! Compute weight
235252
CALL PDAF_local_weight(wtype, rtype, cradius, sradius, distance, 1, 1, tmp, 1.0, weight, 0)
@@ -254,7 +271,23 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH)
254271
dx = abs(clmobs_lon(obs_pdaf2nc(j)) - clmobs_lon(obs_pdaf2nc(i)))
255272
dy = abs(clmobs_lat(obs_pdaf2nc(j)) - clmobs_lat(obs_pdaf2nc(i)))
256273

257-
distance = sqrt(real(dx)**2 + real(dy)**2)
274+
! Check for longitude differences that may yield differences
275+
! larger than 180deg depending on conventions. Example:
276+
! crossing the prime meridian (lon=0deg), when convention is
277+
! all-positive lons
278+
IF (dx > 180.0) THEN
279+
dx = 360.0 - dx
280+
END IF
281+
282+
! Intermediate latitude
283+
yhalf = (clmobs_lat(obs_pdaf2nc(j)) + clmobs_lat(obs_pdaf2nc(i))) / 2.0
284+
285+
! Latitude-dependent factor for longitude difference
286+
dx = dx * cos(yhalf * 3.14159265358979323846 / 180.0)
287+
288+
! Factor ca. 111km comes from R*pi/180, where R is earth
289+
! radius and pi/180 is because we input lat/lon in degrees
290+
distance = 111.19492664455873 * sqrt(real(dx)**2 + real(dy)**2)
258291

259292
! Compute weight
260293
CALL PDAF_local_weight(wtype, rtype, cradius, sradius, distance, 1, 1, tmp, 1.0, weight, 0)
@@ -264,6 +297,15 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH)
264297

265298
END DO
266299
END DO
300+
301+
! Deallocate lon/lat arrays
302+
303+
! For other filters these are deallocated in clean_obs_nc, invoked
304+
! at the end of init_dim_obs_pdaf. For LEnKF, deallocation is
305+
! moved to end of localize_covar_pdaf as these arrays are still
306+
! used here.
307+
if(allocated(clmobs_lon))deallocate(clmobs_lon)
308+
if(allocated(clmobs_lat))deallocate(clmobs_lat)
267309

268310
ENDIF ! model==tag_model_clm
269311
#endif

‎bldsva/intf_DA/pdaf/framework/mod_read_obs.F90‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -506,6 +506,9 @@ end subroutine get_obsindex_currentobsfile
506506
!> This subroutine deallocates the observation arrays used in
507507
!> subroutine `read_obs_nc`.
508508
subroutine clean_obs_nc()
509+
510+
USE mod_assimilation, ONLY: filtertype
511+
509512
implicit none
510513
! if(allocated(idx_obs_nc))deallocate(idx_obs_nc)
511514
if(allocated(pressure_obs))deallocate(pressure_obs)
@@ -514,8 +517,11 @@ subroutine clean_obs_nc()
514517
!if(allocated(y_idx_obs_nc))deallocate(y_idx_obs_nc)
515518
!if(allocated(z_idx_obs_nc))deallocate(z_idx_obs_nc)
516519
!kuw: clean clm observations
517-
if(allocated(clmobs_lon))deallocate(clmobs_lon)
518-
if(allocated(clmobs_lat))deallocate(clmobs_lat)
520+
IF (.NOT. filtertype == 8) THEN
521+
! For LEnKF lat/lon are used
522+
if(allocated(clmobs_lon))deallocate(clmobs_lon)
523+
if(allocated(clmobs_lat))deallocate(clmobs_lat)
524+
END IF
519525
if(allocated(clm_obs))deallocate(clm_obs)
520526
if(allocated(clmobs_layer))deallocate(clmobs_layer)
521527
if(allocated(clmobs_dr))deallocate(clmobs_dr)

‎bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,7 @@ module enkf_clm_mod
7676
real(r8) :: dtime ! time step increment (sec)
7777
integer :: ier ! error code
7878

79-
character(c_char),bind(C,name="outdir"),target :: outdir
79+
character(kind=c_char,len=100),bind(C,name="outdir"),target :: outdir
8080

8181
logical :: log_print ! true=> print diagnostics
8282
real(r8) :: eccf ! earth orbit eccentricity factor

0 commit comments

Comments
 (0)