source: branches/GRISLIv3/SOURCES/lake_rsl_mod.f90

Last change on this file was 481, checked in by aquiquet, 5 months ago

Cleaning branch: unused variables removed following strict compilation
File size: 35.2 KB
Line 
1!*************************************************************************************************************************************
2! This script creates runoff directions from a topographic dataset based on the algorithm of Wang and Liu (2006) described in:
3! An efficient method for identifying and filling surface depressions in digital elevation model for hydrologic analysis and modelling
4! International Journal of Geographical Information Science, Vol. 20, No. 2, Feb. 2006, 193-213
5
6!* It avoids endroheic basins and provide consistent runoff directions for use in the IPSL model.
7!  It compares favourably with directions created with the 'Fill Sinks' algorithm implemented in
8!  GIS softwares (like ArcGIS or QGIS) which is also based on the Wang and Liu algorithm.
9!  There are some differences in runoff directions related to the way I have interpreted some
10!  specific cases that occur. Basins generated by these directions are thus not identical but
11!  major basins are relatively similar.
12!  The corrected topography is however rigourously similar to that obtained from GIS softwares,
13!  at the exception of the Antarctic continent where I noted some inconsistencies in these
14!  softwares due to lack of east-west periodicity and of understanding the increasing longitudinal resolution.
15
16!* INPUTS:
17!  - a topographic dataset
18
19!* OUTPUTS:
20!  - rtm: runoff directions
21!  - topo: initial topographie
22!  - orog_lake : corrected topographic dataset to avoid endorheic basin
23!  - cmsk: coastal mask
24
25!*************************************************************************************************************************************
26
27!* C Dumas - 26.11.2020
28
29module lake_rsl
30
31  use geography, only:nx,ny
32
33  !$ USE OMP_LIB
34
35  implicit none
36
37  integer :: hi_res_runoff ! select 1 : interpolation on hi resolution topo or 0: no interpolation
38
39  !* Routing code parameter
40  integer, dimension(8,2) :: inc
41  integer,dimension(80,2) :: inclake ! 80 voisins (+4 -4) en i,j
42  real, dimension(nx,ny) :: sealevel_2D_loc
43
44
45  !~ ! topo haute resolution (10km) pour calcul des lacs
46  !~ ! nbr de points topo hi res
47  integer, parameter :: nxhi = 964
48  integer, parameter :: nyhi = 964
49  integer, parameter :: ninterpo = 4 ! interpolation avec les 4 points autours
50
51  real, dimension(nxhi,nyhi) :: Bsoc0_hi, S0_hi, H0_hi ! topo haute resolution actuelle (ref)
52  !~ ! anomalie de topo sur la haute resolution
53  !~ !real, dimension(nxhi,nyhi) :: anom_topo_hires
54  !~ ! anomalie de topo
55  real,dimension(nx,ny) :: anom_Bsoc, anom_S, anom_H ! anomalie de topo vs topo ref (low res)
56
57  !~ ! matrice passage 40km vers 10km
58
59  real, dimension(nxhi,nyhi) :: xcc_hi,ycc_hi,xlong_hi,ylat_hi
60  real, dimension(nxhi,nyhi,4) :: poids  ! poids pour passage de low vers hi_res
61  integer,dimension(nxhi,nyhi,4) :: low_2_hi_i, low_2_hi_j ! correspondance grille low vers grille hi_res
62
63  real, dimension(nxhi,nyhi) :: anom_Bsoc_hi, anom_S_hi, anom_H_hi ! anomalie de topo resol std vs resol_hi
64  real, dimension(nxhi,nyhi) :: Bsoc_hi, S_hi, H_hi ! topo haute resolution a l'instant t
65
66
67contains
68
69  !> SUBROUTINE: input_rsl
70  !! Routine permettant d'introduire des variations locales du niveau marin
71  !! par la presence de lacs identifies via un algo de routage sur la topo de GRISLI
72  !>
73  subroutine input_rsl
74
75    use geography, only:dirnameinp
76    use module3D_phy,only:num_param,num_rep_42,xcc,ycc,S,S0,H,H0,Bsoc,Bsoc0
77    use io_netcdf_grisli, only:Read_Ncdf_var
78   
79    implicit none
80
81    namelist/lake_rsl/hi_res_runoff, file_topo_hi_res, coord_topo_hi_res
82
83    !~ !* Debug mode & nc files
84    character(len=100) :: file_topo_hi_res, coord_topo_hi_res
85
86
87    integer :: ii, i, j, k
88    integer :: ios
89    real*8, dimension(:,:), pointer   :: tab2d => null()   !< tableau 2d pointer
90    real, dimension(ninterpo) :: dist_i, dist_j, dist_ij
91
92    rewind(num_param)
93    read(num_param,lake_rsl)
94
95    write(num_rep_42,'(A)')'!___________________________________________________________'
96    write(num_rep_42,'(A)') '&lake_rsl                                    ! nom du bloc '
97    write(num_rep_42,*)
98    write(num_rep_42,*) 'hi_res_runoff = ', hi_res_runoff
99    write(num_rep_42,'(A,A,A)') 'file_topo_hi_res = ', trim(file_topo_hi_res)
100    write(num_rep_42,'(A,A,A)') 'coord_topo_hi_res = ', trim(coord_topo_hi_res)
101    write(num_rep_42,*)'/'
102    write(num_rep_42,*)
103
104
105    !***************************
106    !*** Initialization step ***
107    !***************************
108
109    !~ write(*,*) "Start Initialization step"
110
111    !* The routing code (1er indice : k=1,8, 2eme indice i=1, j=2)
112    inc(1,1) = 0   ! indice i N
113    inc(1,2) = 1   ! indice j N
114    inc(2,1) = 1   ! indice i NE
115    inc(2,2) = 1   ! indice j NE
116    inc(3,1) = 1   ! indice i E
117    inc(3,2) = 0   ! indice j E
118    inc(4,1) = 1   ! indice i SE
119    inc(4,2) = -1  ! indice j SE
120    inc(5,1) = 0   ! indice i S
121    inc(5,2) = -1  ! indice j S
122    inc(6,1) = -1  ! indice i SW
123    inc(6,2) = -1  ! indice j SW
124    inc(7,1) = -1  ! indice i W
125    inc(7,2) = 0   ! indice j W
126    inc(8,1) = -1  ! indice i NW
127    inc(8,2) = 1   ! indice j NW
128
129    !~ print*
130    !~ print*,'rtm directions :'
131    !~ print*,'1 :  N'
132    !~ print*,'2 : NE'
133    !~ print*,'3 :  E'
134    !~ print*,'4 : SE'
135    !~ print*,'5 :  S'
136    !~ print*,'6 : SW'
137    !~ print*,'7 :  W'
138    !~ print*,'8 : NW'
139    !~ print*
140
141
142    !do ii=1,80
143    ii = 1
144    do j = -4,4
145       do i = -4,4
146          if (.not.((i.eq.0).and.(j.eq.0))) then
147             inclake(ii,1)=i
148             inclake(ii,2)=j
149             ii=ii+1
150          endif
151       enddo
152    enddo
153
154
155    if (hi_res_runoff.eq.1) then
156       ! lecture topo hi resolution :
157       call Read_Ncdf_var('Bsoc',trim(dirnameinp)//trim(file_topo_hi_res),tab2d) ! socle
158       Bsoc0_hi(:,:)=tab2d(:,:)
159       call Read_Ncdf_var('S',trim(dirnameinp)//trim(file_topo_hi_res),tab2d)    ! surface
160       S0_hi(:,:)=tab2d(:,:)
161       call Read_Ncdf_var('H',trim(dirnameinp)//trim(file_topo_hi_res),tab2d)    ! epaisseur
162       H0_hi(:,:)=tab2d(:,:)
163
164       ! lecture des coordonnées geographiques
165       ! les coordonnees sont calculees en °dec avec GMT,
166       ! les longitudes sont comprises entre -180 et +180 (negative a l'Ouest de
167       ! Greenwich et positive a l'Est)
168       open(unit=2004,file=trim(dirnameinp)//trim(coord_topo_hi_res),iostat=ios)
169       do k=1,nxhi*nyhi
170          read(2004,*) i,j,xcc_hi(i,j),ycc_hi(i,j),xlong_hi(i,j),ylat_hi(i,j)
171       enddo
172       close(2004)
173
174
175       ! calcul de la matrice de passage de 40km => hi_res (10km)
176       do j=1,nyhi
177          do i=1,nxhi
178             low_2_hi_i(i,j,1) = max(nint((i-1)/4.),1)
179             low_2_hi_i(i,j,2) = min(low_2_hi_i(i,j,1) + 1,nx)
180             low_2_hi_i(i,j,3) = max(nint((i-1)/4.),1)
181             low_2_hi_i(i,j,4) = min(low_2_hi_i(i,j,1) + 1,nx)
182             low_2_hi_j(i,j,1) = max(nint((j-1)/4.),1)
183             low_2_hi_j(i,j,2) = min(low_2_hi_j(i,j,1) + 1,ny)
184             low_2_hi_j(i,j,3) = min(low_2_hi_j(i,j,1) + 1,ny)
185             low_2_hi_j(i,j,4) = max(nint((j-1)/4.),1)
186             ! calcul de la distance entre le pt hi_res et les 4 points low_res autours
187             do k=1,ninterpo
188                dist_i(k)=abs(xcc(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))-xcc_hi(i,j))
189                dist_j(k)=abs(ycc(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))-ycc_hi(i,j))
190                dist_ij(k)=sqrt(dist_i(k)**2+dist_j(k)**2)
191             enddo
192             do k=1,ninterpo
193                poids(i,j,k)=(sum(dist_ij(:))-dist_ij(k))/(3*sum(dist_ij(:)))
194             enddo
195          enddo
196       enddo
197
198       ! calcul de l'anomalie de topo vs la topo de reference :
199       anom_Bsoc(:,:) = Bsoc(:,:) - Bsoc0(:,:)
200       anom_S(:,:) = S(:,:) - S0(:,:)
201       anom_H(:,:) = H(:,:) - H0(:,:)
202
203       ! interpolation 40km => hi_res
204       anom_Bsoc_hi(:,:)=0.
205       anom_S_hi(:,:)=0.
206       anom_H_hi(:,:)=0.
207
208       do k=1,ninterpo
209          do j=1,nyhi
210             do i=1,nxhi
211                anom_Bsoc_hi(i,j) = anom_Bsoc_hi(i,j) + poids(i,j,k)*anom_Bsoc(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
212                anom_S_hi(i,j) = anom_S_hi(i,j) + poids(i,j,k)*anom_S(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
213                anom_H_hi(i,j) = anom_H_hi(i,j) + poids(i,j,k)*anom_H(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
214             enddo
215          enddo
216       enddo
217    endif
218
219  end subroutine input_rsl
220
221
222  !------------------------------------------------------------------------
223  ! subroutine rsl : caclcul du niveau marin local
224  !------------------------------------------------------------------------
225  subroutine rsl
226
227    use module3D_phy,only:BSOC,Bsoc0,S,S0,H,H0,sealevel,sealevel_2d,time,dtmin,ice,mk,debug_3d
228    use param_phy_mod,only : ro,row
229    use netcdf,only: nf90_clobber,nf90_create,nf90_close,nf90_def_dim,nf90_def_var,nf90_put_var,nf90_float,nf90_enddef
230   
231    implicit none
232
233    integer :: i, j, k
234
235    !~ !* Parameters and arrays used in the Initialization step
236    real, dimension(nx,ny) :: orog, orog_lake, lake_level
237    integer, dimension(nx,ny) :: omsk, cmsk
238    integer, dimension(nx,ny) :: cont!, conttrack, contsave, runoff_msk,mask_search, exut
239    integer, dimension(nx,ny) :: ocean
240    integer, dimension(nx,ny) :: pot_lake
241
242    integer :: countcont ! nbr de continents
243    !~ integer :: nbocean, countocean
244
245    !~ !* Parameters and arrays of the runoff calculation
246    integer, dimension(nx,ny) :: rtm
247
248    integer,dimension(nx,ny) :: lake ! 1 si lac, 0 sinon
249    integer, dimension(nx,ny) :: mask_extlake ! extension du lac sous la glace
250    real, dimension(nx,ny) :: orog_extlake ! niveau du lac sur les zones glace a proximité (rayon de 4 pts de grille)
251
252    real, dimension(nx,ny) :: orog_lake_tmp ! pour calculer orog_lake a partir de orog_lake_hi
253
254    real, dimension(nxhi,nyhi) :: orog_hi ! orographie hi_res
255    real, dimension(nxhi,nyhi) :: S_interp_hi, H_interp_hi ! interpolation de S et H low res vers hi_res sans anomalie (pour faire topo englacee non bruitee)
256    integer, dimension(nxhi,nyhi) :: omsk_hi, cmsk_hi
257    integer, dimension(nxhi,nyhi) :: cont_hi ! masque des points continent hi_res
258    integer, dimension(nxhi,nyhi) :: ocean_hi ! masque des points ocean hi_res
259    integer, dimension(nxhi,nyhi) :: pot_lake_hi
260    integer, dimension(nxhi,nyhi) :: lake_hi ! 1 si lac, 0 sinon
261    integer, dimension(nxhi,nyhi) :: rtm_hi ! runoff direction (1 to 8)
262    real, dimension(nxhi,nyhi) :: lake_level_hi ! niveau du lac sur grille hi_res
263    real, dimension(nxhi,nyhi) :: orog_lake_hi  ! orographie a la surface des lacs
264
265    integer :: countcont_hi ! nbr de continents sur grille hi_res
266
267    ! pour debug : ecriture resultats hi_res
268    character(len=100) :: ncfilecont    ! nom fichier netcdf de sortie
269    character(len=1),dimension(2) :: dimname ! nom des dimensions du fichier netcdf de sortie
270    integer :: ncid, xdimid, ydimid, S_hivarid, H_hivarid, Bsoc_hivarid, rtm_hivarid, cont_hivarid, ocean_hivarid, cmsk_hivarid, orog_lake_hivarid, lake_hivarid
271    integer :: status ! erreur operations netcdf
272
273    ncfilecont = "cont-topo_hi_res.nc"
274    dimname(1)='x'
275    dimname(2)='y'
276
277    !~ print*,'Debut rsl time = ', time
278    ! orography sur grille std :
279
280
281    where ((BSOC(:,:)+H(:,:)*ro/row -sealevel).LT.0.)
282       orog(:,:) = BSOC(:,:) ! le point flotte et socle sous le nivx de la mer
283    elsewhere
284       orog(:,:) = BSOC(:,:) + H(:,:)
285    endwhere
286
287    orog_lake(:,:)=orog(:,:) ! orographie a la surface des lacs
288
289    where (orog(:,:).GT.sealevel) ! les depressions continentales seront ajoutees a ce masque plus loin
290       omsk(:,:) = 1
291    elsewhere
292       omsk(:,:) = 0
293    endwhere
294
295    !~ print*,'1/ rsl'
296
297    if (hi_res_runoff.eq.1) then
298       ! calcul de l'anomalie de topo vs la topo de reference :
299       anom_Bsoc(:,:) = Bsoc(:,:) - Bsoc0(:,:)
300       anom_S(:,:) = S(:,:) - S0(:,:)
301       anom_H(:,:) = H(:,:) - H0(:,:)
302
303       anom_Bsoc_hi(:,:) = 0.
304       anom_S_hi(:,:) = 0.
305       S_interp_hi(:,:) = 0.
306       H_interp_hi(:,:) = 0.
307       anom_H_hi(:,:) = 0.
308
309       ! calcul de la topo hi_res :
310       do k=1,ninterpo
311          do j=1,nyhi
312             do i=1,nxhi
313                anom_Bsoc_hi(i,j) = anom_Bsoc_hi(i,j) + poids(i,j,k)*anom_Bsoc(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
314                anom_S_hi(i,j) = anom_S_hi(i,j) + poids(i,j,k)*anom_S(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
315                S_interp_hi(i,j) = S_interp_hi(i,j) + poids(i,j,k)*S(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
316                H_interp_hi(i,j) = H_interp_hi(i,j) + poids(i,j,k)*H(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
317                !        anom_H_hi(i,j) = anom_H_hi(i,j) + poids(i,j,k)*anom_H(low_2_hi_i(i,j,k),low_2_hi_j(i,j,k))
318             enddo
319          enddo
320       enddo
321
322       Bsoc_hi(:,:) = Bsoc0_hi(:,:) + anom_Bsoc_hi(:,:)
323
324       where (H_interp_hi(:,:).GT.1.) ! zones englacees
325          where (Bsoc_hi(:,:).GT.S_interp_hi(:,:)) ! socle qui sort de la glace
326             S_hi(:,:) = Bsoc_hi(:,:)
327             H_hi(:,:) = 0.
328          elsewhere ! socle sous la glace
329             S_hi(:,:) = S_interp_hi(:,:)
330             H_hi(:,:) = S_interp_hi(:,:) - Bsoc_hi(:,:)
331          endwhere
332       elsewhere ! zones non englacees
333          S_hi(:,:) = max(Bsoc_hi(:,:),sealevel)
334          H_hi(:,:) = 0.
335       endwhere
336
337       ! orography sur grille hi_res :
338       where ((Bsoc_hi(:,:)+H_hi(:,:)*ro/row -sealevel).LT.0.)
339          orog_hi(:,:) = Bsoc_hi(:,:) ! le point flotte et socle sous le nivx de la mer
340       elsewhere
341          orog_hi(:,:) = S_hi(:,:)    ! point terre
342       endwhere
343
344       orog_lake_hi(:,:)=orog_hi(:,:) ! orographie a la surface des lacs
345
346       where (orog_hi(:,:).GT.sealevel) ! les depressions continentales seront ajoutees a ce masque plus loin
347          omsk_hi(:,:) = 1
348       elsewhere
349          omsk_hi(:,:) = 0
350       endwhere
351
352    endif
353
354    !~ print*,'2/ rsl'
355    ! mise à jour du niveau des lacs et du routage :
356    if (mod(abs(TIME),500.).lt.dtmin) then
357       call mask_orog_ocean(nx,ny,omsk,cmsk,cont,ocean,countcont,pot_lake)
358
359       !~   print*,'3/ rsl'
360       if (hi_res_runoff.eq.1) then
361          call mask_orog_ocean(nxhi,nyhi,omsk_hi,cmsk_hi,cont_hi,ocean_hi,countcont_hi,pot_lake_hi)
362          call runoff(nxhi,nyhi,countcont_hi,cont_hi,cmsk_hi,omsk_hi,pot_lake_hi,lake_hi,rtm_hi,lake_level_hi,orog_lake_hi)
363          ! interpolation des variables sur grille std :
364          orog_lake(:,:) = 0.
365          orog_lake_tmp(:,:) = 0.
366          lake(:,:) = 0
367          do j=1,ny
368             do i=1,nx
369                orog_lake_tmp(i,j) = max(sum(orog_lake_hi(i*4-3:i*4,j*4-3:j*4)) / 16., orog(i,j))
370                !        lake(i,j) = nint(sum(lake_hi(i*4-3:i*4,j*4-3:j*4)) / 16.)  ! lac si au moins 50% des points sont lacs
371                lake(i,j) = floor(sum(lake_hi(i*4-3:i*4,j*4-3:j*4)) / 16.)  ! lac si 100% des points hi_res sont lacs
372             enddo
373          enddo
374          where (lake(:,:).EQ.1) orog_lake(:,:)=orog_lake_tmp(:,:)
375
376       else
377          call runoff(nx,ny,countcont,cont,cmsk,omsk,pot_lake,lake,rtm,lake_level,orog_lake)
378       endif
379
380       ! pour eviter la prise en compte des lacs dans les depressions a la surface de la calotte :
381       where ((lake(:,:).eq.1).and.(ice.eq.1))
382          lake(:,:) = 0
383          orog_lake(:,:) = orog(:,:)
384       endwhere
385
386
387       !~   print*,'3/ rsl'
388       call extlake(nx,ny,lake,ice,mk,orog_lake,inclake,mask_extlake,orog_extlake)
389
390       !~   print*,'4/ rsl'
391
392       sealevel_2D_loc(:,:) = sealevel_2D(:,:) ! pour debug : A SUPPRIMER QUAND ON REPASSERA SUR SEALEVEL_2D
393       where (lake(:,:).eq.1)
394          sealevel_2D_loc(:,:)=orog_lake(:,:)
395       elsewhere
396          sealevel_2D_loc(:,:)=sealevel
397       endwhere
398
399       where ((mask_extlake.LT.9999).AND.(ice.EQ.1).AND.(mk.EQ.0).AND.(lake.EQ.0).AND.(BSOC.LT.orog_extlake))
400          sealevel_2D_loc(:,:) = orog_extlake(:,:)
401       endwhere
402
403       if (hi_res_runoff.eq.1) then   !afq-- this block is only for debug as it is called each timestep!!
404
405          !* Save mode
406          status = nf90_create(trim(ncfilecont),NF90_CLOBBER, ncid)
407          !  status = nf90_close(ncid)
408          status = nf90_def_dim(ncid, "x", nxhi, xdimid)
409          status = nf90_def_dim(ncid, "y", nyhi, ydimid)
410          status = nf90_def_var(ncid, "S_hi", nf90_float, (/ xdimid, ydimid /), S_hivarid)
411          status = nf90_def_var(ncid, "H_hi", nf90_float, (/ xdimid, ydimid /), H_hivarid)
412          status = nf90_def_var(ncid, "Bsoc_hi", nf90_float, (/ xdimid, ydimid /), Bsoc_hivarid)
413          status = nf90_def_var(ncid, "rtm_hi", nf90_float, (/ xdimid, ydimid /), rtm_hivarid)
414          status = nf90_def_var(ncid, "cont_hi", nf90_float, (/ xdimid, ydimid /), cont_hivarid)
415          status = nf90_def_var(ncid, "ocean_hi", nf90_float, (/ xdimid, ydimid /), ocean_hivarid)
416          status = nf90_def_var(ncid, "cmsk_hi", nf90_float, (/ xdimid, ydimid /), cmsk_hivarid)
417          status = nf90_def_var(ncid, "orog_lake_hi", nf90_float, (/ xdimid, ydimid /), orog_lake_hivarid)
418          status = nf90_def_var(ncid, "lake_hi", nf90_float, (/ xdimid, ydimid /), lake_hivarid)
419
420          status = nf90_enddef(ncid)
421
422          status = nf90_put_var(ncid, S_hivarid, S_hi)
423          status = nf90_put_var(ncid, H_hivarid, H_hi)
424          status = nf90_put_var(ncid, Bsoc_hivarid, Bsoc_hi)
425          status = nf90_put_var(ncid, rtm_hivarid, rtm_hi)
426          status = nf90_put_var(ncid, cont_hivarid, cont_hi)
427          status = nf90_put_var(ncid, ocean_hivarid, ocean_hi)
428          status = nf90_put_var(ncid, cmsk_hivarid, cmsk_hi)
429          status = nf90_put_var(ncid, orog_lake_hivarid, orog_lake_hi)
430          status = nf90_put_var(ncid, lake_hivarid, lake_hi)
431
432          status = nf90_close(ncid)
433       endif
434
435
436    endif
437
438    sealevel_2D(:,:) = sealevel_2D_loc(:,:)
439
440    debug_3d(:,:,126)=sealevel_2D_loc(:,:)
441    debug_3d(:,:,127)=rtm(:,:)
442    debug_3d(:,:,128)=lake(:,:)
443
444    !~ print*,'Fin rsl'
445
446  end subroutine rsl
447
448
449
450  !----------------------------------------------------------------------------
451  ! Subroutine de recherche des points terre, cotiers et ocean
452  !----------------------------------------------------------------------------
453  subroutine mask_orog_ocean(nx_loc,ny_loc,omsk_loc,cmsk_loc,cont_loc,ocean_loc,countcont_loc,pot_lake_loc)
454
455    implicit none
456
457    integer,intent(in) :: nx_loc,ny_loc
458    integer,dimension(nx_loc,ny_loc), intent(inout) :: omsk_loc
459    integer,dimension(nx_loc,ny_loc), intent(out) :: cmsk_loc, cont_loc
460    integer, dimension(nx_loc,ny_loc), intent(out) :: ocean_loc
461    integer, intent(out) :: countcont_loc
462    integer, dimension(nx_loc,ny_loc), intent(out) :: pot_lake_loc
463
464
465    integer :: undefint=0
466    integer, dimension(nx_loc,ny_loc) :: conttrack
467    integer :: sumcont ! somme des cont des 8 points autour du point ii jj
468    integer :: maxloop ! pour debug evite boucle infinies
469    integer, dimension(8) :: ipinctab, jpinctab ! indice des points a scanner autour du point ii jj
470    integer, dimension(8) :: loccont ! cont des 8 points autour du point ii jj
471    integer, dimension(nx_loc,ny_loc) :: contsave
472
473    integer :: nbocean, countocean
474    integer :: countpot_lake ! nbr de dépressions (potential lakes)
475    logical :: stop_lake
476    integer, dimension(2) :: ij_lake
477
478    integer,dimension(8) :: locomsk ! omsk des 8 voisins
479    integer :: maxcont ! valeur max de cont parmi les 8 voisins
480    integer :: sumocean ! somme des ocean des 8 points autour du point ii jj
481    integer, dimension(8) :: lococean ! ocean des 8 points autour du point ii jj
482    integer :: maxocean ! valeur max de ocean des 8 points autour du point ii jj
483    integer,dimension(nx_loc,ny_loc) :: oceansave, true_ocean
484    integer, parameter :: oceanmax=3000 ! nbr max de zones ocean ou lac
485    logical, dimension(oceanmax) :: k_ocean ! permet d'identifier les ocean ouverts des lacs
486    integer :: nbcont
487    integer :: kloc
488    integer :: i,j,ii,jj,k,l
489
490    !* Initialize the coastal mask (cmsk_loc)
491    cmsk_loc = 0
492    !~ print*,'1/ mask_orog_ocean'
493
494    !* Discriminate continents
495    !~ write(*,*) "Number the different continents"
496    cont_loc = undefint
497    nbcont = 1
498    conttrack = 0
499
500    do j = 1, ny_loc
501       do i = 1, nx_loc
502          !do j = 110,111
503          !  do i = 182,183
504
505          !* Locate current grid point
506          ii = i; jj = j
507          sumcont=0
508          maxloop=0
509
510          !* Only on continental points
511          if ( omsk_loc(ii,jj) .eq. 1) then
512             do k=1,8
513                ipinctab(k) = max(1,min(nx_loc,ii+inc(k,1)))
514                ! cas cyclicité selon i :
515                !        ipinctab(k) = ii+inc(k,1)
516                !        if (ipinctab(k).eq.(nx + 1)) ipinctab(k) = 1
517                !        if (ipinctab(k).eq.(0)) ipinctab(k) = nx
518                jpinctab(k) = max(1,min(ny_loc,jj+inc(k,2)))
519                loccont(k)=cont_loc(ipinctab(k),jpinctab(k))
520                sumcont = sumcont + cont_loc(ipinctab(k),jpinctab(k))
521             enddo
522
523             if (sumcont.GT.0) then
524                cont_loc(ii,jj)=minval(loccont,mask=loccont>0)
525                do while (maxval(loccont,mask=loccont>0).GT.cont_loc(ii,jj).and.maxloop.lt.8) ! test si dans les 8 voisins certains ont un loccont > minval
526                   maxcont=maxval(loccont,mask=loccont>0)
527                   where (cont_loc(:,:).EQ.maxcont) cont_loc(:,:) = cont_loc(ii,jj)
528                   where (loccont(:).EQ.maxcont) loccont(:) = cont_loc(ii,jj)
529                   maxloop = maxloop + 1
530                enddo
531             else
532                cont_loc(ii,jj) = nbcont
533                nbcont = nbcont + 1
534             endif
535             conttrack(ii,jj) = 1
536          endif
537       enddo
538    enddo
539    !~ print*,'2/ mask_orog_ocean'   
540
541    !* Reordering to get an incremental list
542    countcont_loc = 1
543    contsave(:,:)=cont_loc(:,:) ! copie cont for reordonation
544
545    if (nbcont.gt.0) then ! on a trouve des continents !!!!
546       do k=1,nbcont
547          if (count(contsave(:,:).eq.k).gt.0) then ! there is continents associated with number k
548             where (cont_loc(:,:).eq.k) cont_loc(:,:) = countcont_loc
549             countcont_loc = countcont_loc + 1
550          endif
551       enddo
552    else
553       print*,'CONTINENT NOT FOUND'
554    endif
555    !~ print*,'3/ mask_orog_ocean'
556
557    ! recherche des lacs et oceans sous le niveau de la mer
558    !~ write(*,*) "Number the different lakes & oceans"
559    ocean_loc = undefint
560    nbocean = 1
561
562    do j = 1, ny_loc
563       do i = 1, nx_loc
564          !* Locate current grid point
565          ii = i; jj = j
566          sumocean=0
567          maxloop=0
568
569          !* Only on ocean points
570          if ( omsk_loc(ii,jj) .eq. 0) then  ! omask=0 => orog < sealevel
571             do k=1,8
572                ipinctab(k) = max(1,min(nx_loc,ii+inc(k,1)))
573                ! cas cyclicité selon i :
574                !        ipinctab(k) = ii+inc(k,1)
575                !        if (ipinctab(k).eq.(nx + 1)) ipinctab(k) = 1
576                !        if (ipinctab(k).eq.(0)) ipinctab(k) = nx
577                jpinctab(k) = max(1,min(ny_loc,jj+inc(k,2)))
578                lococean(k)=ocean_loc(ipinctab(k),jpinctab(k))
579                sumocean = sumocean + ocean_loc(ipinctab(k),jpinctab(k))
580             enddo
581
582             if (sumocean.GT.0) then
583                ocean_loc(ii,jj)=minval(lococean,mask=lococean>0)
584                do while (maxval(lococean,mask=lococean>0).GT.ocean_loc(ii,jj).and.maxloop.lt.8) ! test si dans les 8 voisins certains ont un lococean > minval
585                   maxocean=maxval(lococean,mask=lococean>0)
586                   where (ocean_loc(:,:).EQ.maxocean) ocean_loc(:,:) = ocean_loc(ii,jj)
587                   where (lococean(:).EQ.maxocean) lococean(:) = ocean_loc(ii,jj)
588                   maxloop = maxloop + 1
589                enddo
590             else
591                ocean_loc(ii,jj) = nbocean
592                nbocean = nbocean + 1
593             endif
594          endif
595       enddo
596    enddo
597
598    !~ print*,'4/ mask_orog_ocean' 
599
600    !* Reordering to get an incremental list
601    countocean = 1
602    countpot_lake = 1 ! nombre de dépression intérieur (lacs potentiels
603    oceansave(:,:)=ocean_loc(:,:) ! copie ocean for reordonation
604    true_ocean(:,:)=0 ! point vraiment ocean
605    pot_lake_loc(:,:)=0
606
607    k_ocean(:)=.false.
608
609
610    if (nbocean.gt.0) then ! on a trouve des oceans !!!!
611       ! si l'ocean a un point qui touche un des bords de grille c'est bien un ocean, sinon c'est un lac
612       do i=1,nx_loc
613          if (oceansave(i,1).GT.0) k_ocean(oceansave(i,1))=.true.
614          if (oceansave(i,ny_loc).GT.0) k_ocean(oceansave(i,ny_loc))=.true.
615       enddo
616       do j=1,ny_loc
617          if (oceansave(1,j).GT.0) k_ocean(oceansave(1,j))=.true.
618          if (oceansave(nx_loc,j).GT.0) k_ocean(oceansave(nx_loc,j))=.true.
619       enddo
620
621       do k=1,nbocean
622          if (count(oceansave(:,:).eq.k).gt.0) then ! there is oceans points associated with number k
623             if (k_ocean(k)) then  ! point is true ocean
624                where (oceansave(:,:).eq.k)
625                   ocean_loc(:,:) = countocean
626                endwhere
627                countocean = countocean + 1
628             else ! point is a depresion under sea level
629                where (oceansave(:,:).eq.k)
630                   pot_lake_loc(:,:) = countpot_lake
631                   ocean_loc(:,:) = undefint
632                endwhere
633                countpot_lake = countpot_lake + 1
634             endif
635          endif
636       enddo
637    else
638       print*,'OCEAN NOT FOUND'
639    endif
640    !~ print*,'5/ mask_orog_ocean countpot_lake', countpot_lake
641    ! les points pot_lake sont continent => on leur attribue le no du continent adjacent.
642
643    do l=1,countpot_lake
644       stop_lake = .true.
645       ij_lake=minloc(pot_lake_loc(:,:),mask=(pot_lake_loc.EQ.l))
646       do while (stop_lake)
647          do k=1,8
648             ipinctab(k) = max(1,min(nx_loc,ij_lake(1)+inc(k,1)))
649             ! si cyclicite selon i :
650             !      ipinctab(k) = ijcoord(1)+inc(k,1)
651             !      if (ipinctab(k).eq.(nx + 1)) ipinctab(k) = 1
652             !      if (ipinctab(k).eq.(0)) ipinctab(k) = nx
653             jpinctab(k) = max(1,min(ny_loc,ij_lake(2)+inc(k,2)))
654             locomsk(k)=omsk_loc(ipinctab(k),jpinctab(k))
655             loccont(k)=cont_loc(ipinctab(k),jpinctab(k))
656          enddo
657          if (maxval(loccont).GT.0) then ! voisin point terre
658             ! le lac appartient a ce continent
659             kloc=maxloc(loccont,dim=1)
660             where (pot_lake_loc(:,:).EQ.l)
661                cont_loc(:,:)=cont_loc(ipinctab(kloc),jpinctab(kloc))
662                omsk_loc(:,:)=1
663             endwhere
664             stop_lake=.false.
665          else
666             ij_lake(1) = ij_lake(1)+1
667          endif
668       enddo
669    enddo
670    !~ print*,'6/ mask_orog_ocean'
671
672    ! recherche des points cotiers
673    where ((omsk_loc(:,:).EQ.1).AND.((eoshift(omsk_loc,shift=-1,dim=1).EQ.0) .OR. &
674         (eoshift(omsk_loc,shift=1,dim=1).EQ.0) .OR. (eoshift(omsk_loc,shift=-1,dim=2).EQ.0) .OR. (eoshift(omsk_loc,shift=1,dim=2).EQ.0)))
675       cmsk_loc(:,:) = 1
676    endwhere
677    !~ print*,'7/ mask_orog_ocean'
678
679  end subroutine mask_orog_ocean
680
681
682
683  !---------------------------------------------------------
684  ! subroutine qui calcul l'écoulement (runoff) et les lacs
685  !---------------------------------------------------------
686  subroutine runoff(nx_loc,ny_loc,countcont_loc,cont_loc,cmsk_loc,omsk_loc,pot_lake_loc,lake_loc,rtm_loc,lake_level_loc,orog_lake_loc)
687
688    implicit none
689
690    integer, intent(in) :: nx_loc, ny_loc
691    integer, intent(in) :: countcont_loc ! nbr de continents
692
693    integer,dimension(nx_loc,ny_loc), intent(in) :: cont_loc
694    integer,dimension(nx_loc,ny_loc), intent(in) :: cmsk_loc
695    integer,dimension(nx_loc,ny_loc), intent(in) :: omsk_loc
696    integer,dimension(nx_loc,ny_loc), intent(in) :: pot_lake_loc
697
698    integer,dimension(nx_loc,ny_loc), intent(out) :: lake_loc ! 1 si lac, 0 sinon
699    integer,dimension(nx_loc,ny_loc), intent(out) :: rtm_loc ! runoff direction (1 to 8)
700    real,dimension(nx_loc,ny_loc), intent(out) :: lake_level_loc ! lake level
701    real,dimension(nx_loc,ny_loc), intent(inout) :: orog_lake_loc  ! orographie a la surface des lacs
702
703    integer,dimension(nx_loc,ny_loc) :: exut ! > 0 si exutoire, 0 sinon
704
705    integer,dimension(nx_loc,ny_loc) :: mask_search
706    integer,dimension(nx_loc,ny_loc) :: runoff_msk
707    integer :: nbcont ! no du continent étudié
708    integer,dimension(2) :: ijcoord ! coordonnees ij du point le plus bas
709
710    integer,dimension(2) :: ijcoord_cmsk ! coordonnees ij du point cotier le plus bas
711    integer :: contexut ! compte les points sous le niveau du lac
712    real,dimension(8) :: locorog ! orog des 8 points voisins
713    integer,dimension(8) :: locexut ! 1 si exutoire => permet de gerer les lacs
714    integer,dimension(8) :: ipinctab,jpinctab ! indice des points a scanner autour du point ii jj
715    integer :: countpot_lake ! nbr de dépressions (potential lakes)
716
717
718    integer :: k,l
719
720    !**************************
721    !*** Runoff calculation ***
722    !**************************
723
724    write(*,*) "Start runoff calculation"
725
726
727!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
728    ! version inspiree article
729    !$OMP PARALLEL
730    !$OMP WORKSHARE
731    mask_search(:,:)=0 ! masque des points a prendre en compte a l'iteration n : 0 pas traité, 1 en cours, 2 traité
732    runoff_msk(:,:)=0  ! initialisation du mask points non traités => 0
733    lake_loc(:,:)=0
734    exut(:,:)=0 ! point exutoire d'un lac => 1
735    !output(:,:)=0 ! initialisation des points de sortie vers la mer (0 pas sortie, 1 pt de sortie)
736    !* We work on each continent separately for convenience
737    !do nbcont = 1,1
738    !$OMP END WORKSHARE
739    !$OMP END PARALLEL
740
741    do nbcont = 1,countcont_loc
742       write(*,*) "Working on continent nb :", nbcont
743       ! on demarre des points cotiers en partant du plus bas puis en remontant en altitude
744       ! mise a jour du masque des points a etudier :
745
746       ! demarrage par les points cotier du continent scanné :
747       where ((cont_loc.eq.nbcont).and.(cmsk_loc.eq.1)) mask_search(:,:)=1 ! masque de scan : points cote pour demarrer
748
749       ! debut de la boucle de recherche du point a scanner
750       do while (minval(runoff_msk(:,:),mask=mask_search(:,:).EQ.1).EQ.0) ! test si il reste des points a scanner
751          ijcoord(:)=minloc(orog_lake_loc,mask=(mask_search(:,:).EQ.1))   ! position du point le plus bas etudié
752
753          if (count(mask=((mask_search.eq.1).and.(cmsk_loc.eq.1))).GT.0) then! il reste des points cotiers non traités
754             ! position du point cotier le plus bas :
755             ijcoord_cmsk(:)=minloc(orog_lake_loc,mask=((mask_search(:,:).EQ.1).and.cmsk_loc(:,:).eq.1))
756             ! si altitude egale entre 2 points on traite en priorité le point cotier
757             if (orog_lake_loc(ijcoord(1),ijcoord(2)).EQ.orog_lake_loc(ijcoord_cmsk(1),ijcoord_cmsk(2))) then
758                ijcoord(:) = ijcoord_cmsk(:)
759             endif
760          endif
761          contexut = 0
762          do k=1,8
763             ipinctab(k) = max(1,min(nx_loc,ijcoord(1)+inc(k,1)))
764             ! si cyclicite selon i :
765             !      ipinctab(k) = ijcoord(1)+inc(k,1)
766             !      if (ipinctab(k).eq.(nx + 1)) ipinctab(k) = 1
767             !      if (ipinctab(k).eq.(0)) ipinctab(k) = nx
768             jpinctab(k) = max(1,min(ny_loc,ijcoord(2)+inc(k,2)))
769
770             ! on ajoute les points adjacents aux points scannés pour le prochain tour sauf si ce sont des points ocean ou dejà traités:
771             if ((omsk_loc(ipinctab(k),jpinctab(k)).EQ.1).AND.(mask_search(ipinctab(k),jpinctab(k)).EQ.0)) then
772                mask_search(ipinctab(k),jpinctab(k))=1
773                ! cas des depressions :
774                ! si l'un des points est plus bas que le point de base et n'est pas cotier => on remonte son orog => orog_lake_loc
775                if ((orog_lake_loc(ipinctab(k),jpinctab(k))).LE.(orog_lake_loc(ijcoord(1),ijcoord(2))) &
776                     .AND.(cmsk_loc(ipinctab(k),jpinctab(k)).NE.1)) then
777                   orog_lake_loc(ipinctab(k),jpinctab(k)) = orog_lake_loc(ijcoord(1),ijcoord(2))
778                   lake_loc(ipinctab(k),jpinctab(k)) = 1
779                   contexut= contexut + 1
780                endif
781             endif
782             locorog(k) = orog_lake_loc(ipinctab(k),jpinctab(k))
783             locexut(k) = exut(ipinctab(k),jpinctab(k))
784          enddo
785          if ((contexut.ge.1).and.(exut(ijcoord(1),ijcoord(2)).EQ.0)) exut(ijcoord(1),ijcoord(2)) = 1
786          rtm_loc(ijcoord(1),ijcoord(2))=minloc(locorog,dim=1) ! sens ecoulement donne par minloc de la topo des 8 voisins
787          ! cas d'un point lac
788          if (lake_loc(ijcoord(1),ijcoord(2)).eq.1) then
789             exut(ijcoord(1),ijcoord(2)) = maxval(locexut) + 1
790             ! ecoulement vers l'exutoire du lac
791             if (maxval(locexut).GT.0) then
792                rtm_loc(ijcoord(1),ijcoord(2))=minloc(locexut,dim=1,mask=(locexut.ge.1)) ! ecoulement vers l'exutoire
793             endif
794          endif
795          ! on supprime le point attribué de mask_search :
796          mask_search(ijcoord(1),ijcoord(2))=2
797          ! le point a ete traité :
798          runoff_msk(ijcoord(1),ijcoord(2))=1
799       enddo
800    enddo
801
802    ! pour chaque lac calcul du niveau max :
803    lake_level_loc(:,:)=-9999. ! initialisation
804    do l=1,countpot_lake
805       where (pot_lake_loc(:,:).EQ.l)
806          lake_level_loc(:,:)=orog_lake_loc(:,:)
807       endwhere
808    enddo
809
810  end subroutine runoff
811
812  !--------------------------------------------------------------------
813  !  subroutine extlake : extension des lacs sous la glace
814  !--------------------------------------------------------------------
815
816  subroutine extlake(nx_loc,ny_loc,lake_loc,ice_loc,mk_loc,orog_lake_loc,inclake_loc,mask_extlake_loc,orog_extlake_loc)
817
818    implicit none
819
820    integer, intent(in) :: nx_loc, ny_loc
821    integer,dimension(nx_loc,ny_loc), intent(in) :: lake_loc ! 1 si lac, 0 sinon
822    integer,dimension(nx_loc,ny_loc), intent(in) :: ice_loc ! 1 si glace, 0 sinon
823    integer,dimension(nx_loc,ny_loc), intent(in) :: mk_loc ! 0 si terre, 1 sinon
824    real,dimension(nx_loc,ny_loc), intent(in) :: orog_lake_loc ! orographie a la surface des lacs
825    integer,dimension(80,2),intent(in) :: inclake_loc ! 80 voisins (+4 -4) en i,j
826    integer, dimension(nx_loc,ny_loc),intent(out) :: mask_extlake_loc ! extension du lac sous la glace
827    real, dimension(nx_loc,ny_loc),intent(out) :: orog_extlake_loc ! niveau du lac sur les zones glace a proximité (rayon de 4 pts de grille)
828
829    integer, dimension(8) :: ipinctab,jpinctab ! indice des points a scanner autour du point ii jj
830    integer, dimension(nx_loc,ny_loc) :: clake ! point lac avec voisin point englace
831    integer :: i,j,k
832
833    ! initialisation
834    clake(:,:)=0
835    orog_extlake_loc(:,:)=0.
836    mask_extlake_loc(:,:)=9999
837
838    ! recherche des points lacs en bordure de calotte => voisin pas lac englace et terre
839    where ((lake_loc(:,:).EQ.1).AND.(((eoshift(lake_loc,shift=-1,dim=1).EQ.0).AND.(eoshift(ice_loc,shift=-1,dim=1).EQ.1).AND.(eoshift(mk_loc,shift=-1,dim=1).EQ.0)) .OR. &
840         ((eoshift(lake_loc,shift=1,dim=1).EQ.0).AND.(eoshift(ice_loc,shift=1,dim=1).EQ.1).AND.(eoshift(mk_loc,shift=1,dim=1).EQ.0)) .OR. &
841         ((eoshift(lake_loc,shift=-1,dim=2).EQ.0).AND.(eoshift(ice_loc,shift=-1,dim=2).EQ.1).AND.(eoshift(mk_loc,shift=-1,dim=2).EQ.0)) .OR. &
842         ((eoshift(lake_loc,shift=1,dim=2).EQ.0).AND.(eoshift(ice_loc,shift=1,dim=2).EQ.1).AND.(eoshift(mk_loc,shift=1,dim=2).EQ.0))))
843       clake(:,:) = 1
844    endwhere
845
846
847
848    ! point terre et englace : mk=0 et H>0
849    do j=1,ny
850       do i=1,nx
851          if (clake(i,j).EQ.1) then
852             do k=1,80
853                ipinctab(k) = max(1,min(nx,i + inclake_loc(k,1)))
854                jpinctab(k) = max(1,min(ny,j + inclake_loc(k,2)))
855                mask_extlake_loc(ipinctab(k),jpinctab(k)) = 1
856                orog_extlake_loc(ipinctab(k),jpinctab(k))=orog_lake_loc(i,j)
857             enddo
858          endif
859       enddo
860    enddo
861
862  end subroutine extlake
863
864end module lake_rsl
865
866
867
868
869
870
871
872
Note: See TracBrowser for help on using the repository browser.