Changeset 329

Timestamp:

01/26/21 11:47:38 (3 years ago)

Author:

dumas

Message:

hi-resolution lake update in lake_rsl_mod.f90

File:

• trunk/SOURCES/lake_rsl_mod.f90 (modified) (26 diffs)

trunk/SOURCES/lake_rsl_mod.f90

@@ -172 +172 @@
 !~ write(*,*) "Start Initialization step"
 
-!* The routing code (i=1, j=2)
-inc(1,1) = 0
-inc(1,2) = 1
-inc(2,1) = 1
-inc(2,2) = 1
-inc(3,1) = 1
-inc(3,2) = 0
-inc(4,1) = 1
-inc(4,2) = -1
-inc(5,1) = 0
-inc(5,2) = -1
-inc(6,1) = -1
-inc(6,2) = -1
-inc(7,1) = -1
-inc(7,2) = 0
-inc(8,1) = -1
-inc(8,2) = 1
+!* The routing code (1er indice : k=1,8, 2eme indice i=1, j=2)
+inc(1,1) = 0   ! indice i N
+inc(1,2) = 1   ! indice j N
+inc(2,1) = 1   ! indice i NE
+inc(2,2) = 1   ! indice j NE
+inc(3,1) = 1   ! indice i E
+inc(3,2) = 0   ! indice j E
+inc(4,1) = 1   ! indice i SE
+inc(4,2) = -1  ! indice j SE
+inc(5,1) = 0   ! indice i S
+inc(5,2) = -1  ! indice j S
+inc(6,1) = -1  ! indice i SW
+inc(6,2) = -1  ! indice j SW
+inc(7,1) = -1  ! indice i W
+inc(7,2) = 0   ! indice j W
+inc(8,1) = -1  ! indice i NW
+inc(8,2) = 1   ! indice j NW
+
+!~ print*
+!~ print*,'rtm directions :'
+!~ print*,'1 :  N'
+!~ print*,'2 : NE'
+!~ print*,'3 :  E'
+!~ print*,'4 : SE'
+!~ print*,'5 :  S'
+!~ print*,'6 : SW'
+!~ print*,'7 :  W'
+!~ print*,'8 : NW'
+!~ print*
 
 
@@ -230 +242 @@
 
 ! calcul de la matrice de passage de 40km => hi_res (10km)
-!do j=5,ny-4
-!  do i=5,nx-4
   do j=1,nyhi
     do i=1,nxhi
-!j=5
-!i=5
-!    do k=1,ninterpo
-!      low_2_hi(i,j,k) = nint((i-1)/4)
-!    enddo
       low_2_hi_i(i,j,1) = max(nint((i-1)/4.),1)
       low_2_hi_i(i,j,2) = min(low_2_hi_i(i,j,1) + 1,nx)
@@ -280 +285 @@
 endif
 
-!~ print*,'Input RSL'
-!~ print*,'65,166 S0 H0 Bsoc0 ',S0(65,166), H0(65,166), Bsoc0(65,166)
-!~ print*,'260,664 S0 H0 Bsoc0 ',S0_hi(260,664), H0_hi(260,664), Bsoc0_hi(260,664)
-!~ print*,'65,166 S H Bsoc ',S(65,166), H(65,166), Bsoc(65,166)
-!~ print*,'260,664 Shi Hhi Bsochi ',S_hi(260,664), H_hi(260,664), Bsoc_hi(260,664)
-!~ print*,'260,664 anom Shi Hhi Bsochi ',anom_S_hi(260,664), anom_H_hi(260,664), anom_Bsoc_hi(260,664)
-
-
-!print*, 'input_rsl poids', poids(5,5,:)
-!print*,low_2_hi_i(i,j,:)
-!print*,low_2_hi_j(i,j,:)
-!print*,dist_i(:)
-!print*,dist_j(:)
-!print*,dist_ij(:)
-!print*,sum(dist_ij(:))
-!print*,xcc(1,1),xcc(2,2),ycc(1,1),ycc(2,2),xcc(nx,ny),ycc(nx,ny)
-!print*,xcc_hi(1,1),ycc_hi(1,1)
-!print*,xcc_hi(5,5),xcc_hi(6,6),ycc_hi(5,5),ycc_hi(6,6),xcc_hi(nxhi,nyhi),ycc_hi(nxhi,nyhi)
-!print*, xcc(low_2_hi_i(i,j,1),low_2_hi_j(i,j,1))-xcc_hi(i,j)
-!~ print*,'Fin input_rsl'
-
 end subroutine input_rsl
 
@@ -318 +302 @@
 character(len=100) :: ncfilecont, ncfilernff, ncfiletopoout,ncfilecmsk    ! nom fichier netcdf de sortie
 character(len=1),dimension(2) :: dimname ! nom des dimensions du fichier netcdf de sortie
-integer :: ncid, xdimid, ydimid, S_hivarid, H_hivarid, Bsoc_hivarid, rtm_hivarid, cont_hivarid, ocean_hivarid, cmsk_hivarid, orog_lake_hivarid
+integer :: ncid, xdimid, ydimid, S_hivarid, H_hivarid, Bsoc_hivarid, rtm_hivarid, cont_hivarid, ocean_hivarid, cmsk_hivarid, orog_lake_hivarid, lake_hivarid
 integer :: status ! erreur operations netcdf
 
@@ -370 +354 @@
   enddo
   
-  Bsoc_hi(:,:) = Bsoc_hi(:,:) + anom_Bsoc_hi(:,:)
-  where ((Bsoc_hi(:,:).GT.S_interp_hi(:,:)).and.(H_interp_hi(:,:).GT.1.))
-    S_hi(:,:) = Bsoc_hi(:,:)
+  Bsoc_hi(:,:) = Bsoc0_hi(:,:) + anom_Bsoc_hi(:,:)
+  
+  where (H_interp_hi(:,:).GT.1.) ! zones englacees
+    where (Bsoc_hi(:,:).GT.S_interp_hi(:,:)) ! socle qui sort de la glace
+      S_hi(:,:) = Bsoc_hi(:,:)
+      H_hi(:,:) = 0.
+    elsewhere ! socle sous la glace
+      S_hi(:,:) = S_interp_hi(:,:)
+      H_hi(:,:) = S_interp_hi(:,:) - Bsoc_hi(:,:)
+    endwhere
+  elsewhere ! zones non englacees
+    S_hi(:,:) = max(Bsoc_hi(:,:),sealevel)
     H_hi(:,:) = 0.
-  elsewhere ((Bsoc_hi(:,:) + H_interp_hi(:,:)).GT.S_interp_hi(:,:))
-    S_hi(:,:) = S_interp_hi(:,:)
-    H_hi(:,:) = S_interp_hi(:,:) - Bsoc_hi(:,:)
-  elsewhere (H_interp_hi(:,:).GT.1.)
-    S_hi(:,:) = S_interp_hi(:,:)
-    H_hi(:,:) = max(0., S_interp_hi(:,:) - Bsoc_hi(:,:))
   endwhere
-  
-  
-  
-!~   where ((H_interp_hi(:,:).GT.10.).and.(Bsoc_hi(:,:)+) ! zone englacee => topo S_interp_hi (pas en anomalie pour eviter topo bruitee)
-!~     S_hi(:,:) = S_interp_hi(:,:)
-!~     H_hi(:,:) = max(0., S_interp_hi(:,:) - Bsoc_hi(:,:))
-!~   elsewhere
-!~     S_hi(:,:) = S_hi(:,:) + anom_S_hi(:,:)
-!~     H_hi(:,:) = 0.
-!~   endwhere
-
   
 ! orography sur grille hi_res :
@@ -407 +383 @@
     omsk_hi(:,:) = 0
   endwhere
-
-!~ print*,'RSL Hi_res'
-!~ print*,'65,166 S0 H0 Bsoc0 ',S0(65,166), H0(65,166), Bsoc0(65,166)
-!~ print*,'260,664 S0 H0 Bsoc0 ',S0_hi(260,664), H0_hi(260,664), Bsoc0_hi(260,664)
-!~ print*,'65,166 S H Bsoc ',S(65,166), H(65,166), Bsoc(65,166)
-!~ print*,'260,664 Shi Hhi Bsochi ',S_hi(260,664), H_hi(260,664), Bsoc_hi(260,664)
-!~ print*,'260,664 anom Shi Hhi Bsochi ',anom_S_hi(260,664), anom_H_hi(260,664), anom_Bsoc_hi(260,664)
-!~ print*,'low_2_hi_i', low_2_hi_i(260,664,:)
-!~ print*,'low_2_hi_i', low_2_hi_j(260,664,:)
-!~ print*, 'poids',poids(260,664,:)
-!~ print*,'low_2_hi_i', low_2_hi_i(1,1,:)
-!~ print*,'low_2_hi_i', low_2_hi_j(1,1,:)
-!~ print*, 'poids',poids(1,1,:)
-!~ print*,'low_2_hi_i', low_2_hi_i(964,964,:)
-!~ print*,'low_2_hi_i', low_2_hi_j(964,964,:)
-!~ print*, 'poids',poids(964,964,:)
 
 endif
@@ -458 +418 @@
 !~   print*,'4/ rsl'
 
-!print*, '39,63 pot_lake',lake_level(39,63),orog_lake(39,63),pot_lake(39,63)
-
-
-!~ print*
-!~ print*,'rtm directions :'
-!~ print*,'1 :  N'
-!~ print*,'2 : NE'
-!~ print*,'3 :  E'
-!~ print*,'4 : SE'
-!~ print*,'5 :  S'
-!~ print*,'6 : SW'
-!~ print*,'7 :  W'
-!~ print*,'8 : NW'
-!~ print*
+  sealevel_2D_loc(:,:) = sealevel_2D(:,:) ! pour debug : A SUPPRIMER QUAND ON REPASSERA SUR SEALEVEL_2D
+  where (lake(:,:).eq.1)
+    sealevel_2D_loc(:,:)=orog_lake(:,:)
+  elsewhere
+    sealevel_2D_loc(:,:)=sealevel
+  endwhere
+
+  where ((mask_extlake.LT.9999).AND.(ice.EQ.1).AND.(mk.EQ.0).AND.(lake.EQ.0).AND.(BSOC.LT.orog_extlake))
+    sealevel_2D_loc(:,:) = orog_extlake(:,:)
+  endwhere
+
   if (hi_res_runoff.eq.1) then
-!~   print*,'AVANT SAVE NETCDF'
-!~   print*,'S_hi',S_hi(1,1),S_hi(65,166)
 
 !* Save mode
@@ -489 +444 @@
     status = nf90_def_var(ncid, "cmsk_hi", nf90_float, (/ xdimid, ydimid /), cmsk_hivarid)
     status = nf90_def_var(ncid, "orog_lake_hi", nf90_float, (/ xdimid, ydimid /), orog_lake_hivarid)
+    status = nf90_def_var(ncid, "lake_hi", nf90_float, (/ xdimid, ydimid /), lake_hivarid)
   
     status = nf90_enddef(ncid)
@@ -500 +456 @@
     status = nf90_put_var(ncid, cmsk_hivarid, cmsk_hi)
     status = nf90_put_var(ncid, orog_lake_hivarid, orog_lake_hi)
+    status = nf90_put_var(ncid, lake_hivarid, lake_hi)
   
     status = nf90_close(ncid)
   endif
-!  call write_ncdf_dim('x',trim(ncfilecont),nxhi)
-!  call write_ncdf_dim('y',trim(ncfilecont),nyhi)
-!  tab(:,:)=S_hi(:,:)
-!  call Write_Ncdf_var('S_hi',dimname,TRIM(ncfilecont),tab,1,0,'double')
-!~   tab(:,:)=H_hi(:,:)
-!~   call Write_Ncdf_var('H_hi',dimname,TRIM(ncfilecont),tab,'double')
-!~   tab(:,:)=Bsoc_hi(:,:)
-!~   call Write_Ncdf_var('Bsoc_hi',dimname,TRIM(ncfilecont),tab,'double')
-!~   tab(:,:)=rtm(:,:)
-!~   call Write_Ncdf_var('rtm',dimname,TRIM(ncfilecont),tab,'double')
-!~   tab(:,:)=orog(:,:)
-!~   call Write_Ncdf_var('cont',dimname,TRIM(ncfilecont),tab,'double')
-!~   tab(:,:)=cmsk(:,:)
-!~   call Write_Ncdf_var('cmsk',dimname,TRIM(ncfilecont),tab,'double')
-!~   tab(:,:)=orog_lake(:,:)
-!~   call Write_Ncdf_var('orog_lake',dimname,TRIM(ncfilecont),tab,'double')
-
-!~   status = nf90_create(trim(ncfilecont),and(nf90_write,nf90_64bit_offset),ncid)
-!~   status = nf90_close(ncid)
-!~   call write_ncdf_dim('x',trim(ncfilecont),nx)
-!~   call write_ncdf_dim('y',trim(ncfilecont),ny)
-!~   tab(:,:)=S(:,:)
-!~   call Write_Ncdf_var('S',dimname,TRIM(ncfilecont),tab,'double')
-
-
-!~   print*, '50,138 lake,sealevel_2D,ice,mk,clake,mask_extlake,orog_extlake '
-!~   print*,S(50,138),H(50,138),BSOC(50,138)
-!~   print*,lake(50,138),sealevel_2D(50,138),ice(50,138),mk(50,138),mask_extlake(50,138),orog_extlake(50,138)
-!~   print*,'65,166'
-!~   print*,S(65,166),H(65,166),BSOC(65,166)
-!~   print*,lake(65,166),sealevel_2D(65,166),ice(65,166),mk(65,166),mask_extlake(65,166),orog_extlake(65,166)
-
-
-!~ print*,'153,41'
-!~ print*,S(153,41),H(153,41),BSOC(153,41)
-!~ print*,lake(153,41),sealevel_2D(153,41),ice(153,41),mk(153,41),clake(153,41),mask_extlake(153,41),orog_extlake(153,41)
-
-!where (orog_lake(:,:).gt.orog(:,:))  ! pose pb lorsque topo evolue => test devient faux
-  sealevel_2D_loc(:,:) = sealevel_2D(:,:) ! pour debug : A SUPPRIMER QUAND ON REPASSERA SUR SEALEVEL_2D
-  where (lake(:,:).eq.1)
-    sealevel_2D_loc(:,:)=orog_lake(:,:)
-  elsewhere
-    sealevel_2D_loc(:,:)=sealevel
-  endwhere
-
-!where ((mask_extlake.EQ.1).AND.(ice.EQ.1).AND.(mk.EQ.0).AND.(lake.EQ.0).AND.(BSOC.LT.orog_extlake))
-  where ((mask_extlake.LT.9999).AND.(ice.EQ.1).AND.(mk.EQ.0).AND.(lake.EQ.0).AND.(BSOC.LT.orog_extlake))
-    sealevel_2D_loc(:,:) = orog_extlake(:,:)
-!  sealevel_2D(:,:) = sealevel  + (orog_extlake(:,:)-sealevel) / mask_extlake(:,:) 
-  endwhere
+
 
 endif
 
 sealevel_2D(:,:) = sealevel_2D_loc(:,:)
-
-
-
-!~ print*,'56,125',S(56,125),H(56,125),BSOC(56,125),MK(56,125),ice(56,125),omsk(56,125),orog(56,125),orog_lake(56,125)
 
 debug_3d(:,:,126)=sealevel_2D_loc(:,:)
 debug_3d(:,:,127)=cont(:,:)
 debug_3d(:,:,128)=lake(:,:)
-
-!print*,'lake_rsl   print*,'Debut rsl time = ', time59,119', time, sealevel_2D_loc(59,119),orog_lake(59,119),orog(59,119),sealevel,orog_extlake(59,119)
- 
-!print*,'lake_rsl 6,112'
-!print('(6(f8.1,1x),5(i4,1x),1(f8.1))'),time, sealevel_2D_loc(6,112),orog_lake(6,112),orog(6,112),sealevel,orog_extlake(6,112),cont(6,112),lake(6,112),mask_extlake(6,112),ice(6,112),mk(6,112),bsoc(6,112)
 
 !~ print*,'Fin rsl'
@@ -621 +520 @@
 cmsk_loc = 0
 !~ print*,'1/ mask_orog_ocean'
-!~ print*,'mask_orog_ocean inc = ', inc(1,1), inc(1,2)
-!~ print*,'RSL : update lakes and runoff routing', time
-
-!print*,'2/ Point (182,110) orog, omsk, cmsk, cont', orog(182,110), omsk(182,110), cmsk(182,110), cont(182,110)
+
 !* Discriminate continents
 !~ write(*,*) "Number the different continents"
@@ -653 +549 @@
         sumcont = sumcont + cont_loc(ipinctab(k),jpinctab(k))
       enddo
-      !print*,'ii jj loccont',ii,jj,loccont
-      !print*,'sumcont',sumcont  
+ 
       if (sumcont.GT.0) then
         cont_loc(ii,jj)=minval(loccont,mask=loccont>0)
         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
           maxcont=maxval(loccont,mask=loccont>0)
-!          print*,'ii,jj maxval,cont,maxcont,minval', ii,jj, maxval(loccont,mask=loccont>0), cont(ii,jj),maxcont
-!          write(6,*),'loccont',loccont
-!          where (cont(:,:).EQ.maxval(loccont,mask=loccont>0))
           where (cont_loc(:,:).EQ.maxcont) cont_loc(:,:) = cont_loc(ii,jj)
           where (loccont(:).EQ.maxcont) loccont(:) = cont_loc(ii,jj)
-!          print*,'cont voisin(358,10) :',cont(358,10)
           maxloop = maxloop + 1
         enddo
@@ -675 +566 @@
   enddo
 enddo
-!~ print*,'2/ mask_orog_ocean'
-!print*,'3/ Point (182,110) orog, omsk, cmsk, cont', orog(182,110), omsk(182,110), cmsk(182,110), cont(182,110)
-
-
-  
-     
+!~ print*,'2/ mask_orog_ocean'   
 
 !* Reordering to get an incremental list
@@ -705 +591 @@
 do j = 1, ny_loc
   do i = 1, nx_loc
-!do j = 110,111
-!  do i = 182,183 
-
     !* Locate current grid point
     ii = i; jj = j
@@ -725 +608 @@
         sumocean = sumocean + ocean_loc(ipinctab(k),jpinctab(k))
       enddo
-      !print*,'ii jj loccont',ii,jj,loccont
-      !print*,'sumcont',sumcont  
+
       if (sumocean.GT.0) then
         ocean_loc(ii,jj)=minval(lococean,mask=lococean>0)
         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
           maxocean=maxval(lococean,mask=lococean>0)
-!          print*,'ii,jj maxval,cont,maxcont,minval', ii,jj, maxval(loccont,mask=loccont>0), cont(ii,jj),maxcont
-!          write(6,*),'loccont',loccont
-!          where (cont(:,:).EQ.maxval(loccont,mask=loccont>0))
           where (ocean_loc(:,:).EQ.maxocean) ocean_loc(:,:) = ocean_loc(ii,jj)
           where (lococean(:).EQ.maxocean) lococean(:) = ocean_loc(ii,jj)
-!          print*,'cont voisin(358,10) :',cont(358,10)
           maxloop = maxloop + 1
         enddo
@@ -795 +673 @@
   stop_lake = .true.
   ij_lake=minloc(pot_lake_loc(:,:),mask=(pot_lake_loc.EQ.l))
-!  print*,'lake, i, j : ',l, ij_lake(:)
   do while (stop_lake)
     do k=1,8
@@ -808 +685 @@
     enddo
     if (maxval(loccont).GT.0) then ! voisin point terre
-!      print*,'maxval(loccont)',maxval(loccont)
       ! le lac appartient a ce continent
       kloc=maxloc(loccont,dim=1)
@@ -815 +691 @@
         omsk_loc(:,:)=1
       endwhere
-!      print*,'lake => cont ',l, cont(ij_lake(1),ij_lake(2))
       stop_lake=.false.
     else
@@ -906 +781 @@
   do while (minval(runoff_msk(:,:),mask=mask_search(:,:).EQ.1).EQ.0) ! test si il reste des points a scanner
     ijcoord(:)=minloc(orog_lake_loc,mask=(mask_search(:,:).EQ.1))   ! position du point le plus bas etudié
-!debug
-!    if (ijcoord(1).EQ.286.and.ijcoord(2).eq.26) then
-!      print*,'POINT 286 26',  ijcoord, orog_lake(ijcoord(1),ijcoord(2))
-!      do j=1,ny
-!        do i=1,nx
-!          if (mask_search(i,j).EQ.1) then
-!            print*,'POINT 286 26 mask_search', i,j,mask_search(i,j),orog_lake(i,j)
-!          endif
-!        enddo
-!      enddo
-!    endif
-!fin debug
     
     if (count(mask=((mask_search.eq.1).and.(cmsk_loc.eq.1))).GT.0) then! il reste des points cotiers non traités
@@ -924 +787 @@
       ! si altitude egale entre 2 points on traite en priorité le point cotier
       if (orog_lake_loc(ijcoord(1),ijcoord(2)).EQ.orog_lake_loc(ijcoord_cmsk(1),ijcoord_cmsk(2))) then
-!        print*,'point cotier prio : ',ijcoord_cmsk, orog_lake(ijcoord_cmsk(1),ijcoord_cmsk(2))
-!        print*,'point ijcoord remplace : ',ijcoord, orog_lake(ijcoord(1),ijcoord(2))
         ijcoord(:) = ijcoord_cmsk(:)
-!        pause
       endif
     endif
-!    print*,'ijcoord, orog ',ijcoord(:), orog(ijcoord(1),ijcoord(2))
-    ! orog des 8 points voisins => locorog
     contexut = 0
     do k=1,8
@@ -944 +802 @@
       if ((omsk_loc(ipinctab(k),jpinctab(k)).EQ.1).AND.(mask_search(ipinctab(k),jpinctab(k)).EQ.0)) then
         mask_search(ipinctab(k),jpinctab(k))=1
-!      mask_search(ipinctab(k),jpinctab(k))=min(omsk(ipinctab(k),jpinctab(k)),1) ! => si on est sur un point ocean omsk=0
 ! cas des depressions :
 ! 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
@@ -952 +809 @@
           lake_loc(ipinctab(k),jpinctab(k)) = 1
           contexut= contexut + 1
-!          print*, '!!!!! LAKE !!!!! ', orog(ipinctab(k),jpinctab(k)), orog_lake_loc(ipinctab(k),jpinctab(k))
         endif
       endif
       locorog(k) = orog_lake_loc(ipinctab(k),jpinctab(k))
-!      locmsk(k) = runoff_msk(ipinctab(k),jpinctab(k))
-!      loclake(k) = lake(ipinctab(k),jpinctab(k))
-!      locmsk_search = mask_search(ipinctab(k),jpinctab(k))
       locexut(k) = exut(ipinctab(k),jpinctab(k))
     enddo
@@ -966 +819 @@
     if (lake_loc(ijcoord(1),ijcoord(2)).eq.1) then
       exut(ijcoord(1),ijcoord(2)) = maxval(locexut) + 1
-!      print*,'1/ lake ij', ijcoord(1),ijcoord(2), orog_lake_loc(ijcoord(1),ijcoord(2))
-!      print*,'2/ locorog ',locorog
-!      print*,'3/ locexut ',locexut
     ! ecoulement vers l'exutoire du lac
       if (maxval(locexut).GT.0) then
         rtm_loc(ijcoord(1),ijcoord(2))=minloc(locexut,dim=1,mask=(locexut.ge.1)) ! ecoulement vers l'exutoire
-    
-!      if (maxval(locexut).EQ.1) then
-!        rtm(ijcoord(1),ijcoord(2))=maxloc(locexut,dim=1) ! ecoulement vers l'exutoire
-!      else ! pas de voisin exutoire
-      ! ecoulement idem que voisin lac :
-!        rtm(ijcoord(1),ijcoord(2))=maxloc(locmsk,mask=loclake.EQ.1,dim=1) ! localisation du voisin deja traite (locmsk=1) qui est lac (loclake=1)
       endif
-!      print*,'4/ rtm ',rtm(ijcoord(1),ijcoord(2))
-!      print*,'5/ locmsk ',locmsk
-!    else
-!      print*,'point std traite: ',ijcoord, orog_lake_loc(ijcoord(1),ijcoord(2))
-!      print*,'rtm ', rtm(ijcoord(1),ijcoord(2))
-!      print*,'locorog ', locorog
-!      print*,'cmsk ', cmsk(ijcoord(1),ijcoord(2))
-!      print*,'lake ', lake(ijcoord(1),ijcoord(2))
-!      print*
-      
     endif
     ! on supprime le point attribué de mask_search :
@@ -994 +828 @@
     ! le point a ete traité :
     runoff_msk(ijcoord(1),ijcoord(2))=1
-!    minloc(locorog,mask=((locorog.gt.0).and.(locmsk.eq.0)))
-!    print*, rtm(ijcoord(1),ijcoord(2))
-!    contprint = contprint + 1
-!    if (contprint.eq.20) goto 2020
   enddo
 enddo
@@ -1004 +834 @@
 lake_level_loc(:,:)=-9999. ! initialisation
 do l=1,countpot_lake
-!print*, 'pot_lake :',l
   where (pot_lake_loc(:,:).EQ.l)
     lake_level_loc(:,:)=orog_lake_loc(:,:)
@@ -1010 +839 @@
 enddo
 
-
-
 end subroutine runoff
 
@@ -1020 +847 @@
 subroutine extlake(nx_loc,ny_loc,lake_loc,ice_loc,mk_loc,orog_lake_loc,inclake_loc,mask_extlake_loc,orog_extlake_loc)
 
-implicit none
-
-integer, intent(in) :: nx_loc, ny_loc
-integer,dimension(nx_loc,ny_loc), intent(in) :: lake_loc ! 1 si lac, 0 sinon
-integer,dimension(nx_loc,ny_loc), intent(in) :: ice_loc ! 1 si glace, 0 sinon
-integer,dimension(nx_loc,ny_loc), intent(in) :: mk_loc ! 0 si terre, 1 sinon
-real,dimension(nx_loc,ny_loc), intent(in) :: orog_lake_loc ! orographie a la surface des lacs
-integer,dimension(80,2),intent(in) :: inclake_loc ! 80 voisins (+4 -4) en i,j
-integer, dimension(nx_loc,ny_loc),intent(out) :: mask_extlake_loc ! extension du lac sous la glace
-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)
-
-integer, dimension(8) :: ipinctab,jpinctab ! indice des points a scanner autour du point ii jj
-integer, dimension(nx_loc,ny_loc) :: clake ! point lac avec voisin point englace
-integer :: i,j,k
+  implicit none
+
+  integer, intent(in) :: nx_loc, ny_loc
+  integer,dimension(nx_loc,ny_loc), intent(in) :: lake_loc ! 1 si lac, 0 sinon
+  integer,dimension(nx_loc,ny_loc), intent(in) :: ice_loc ! 1 si glace, 0 sinon
+  integer,dimension(nx_loc,ny_loc), intent(in) :: mk_loc ! 0 si terre, 1 sinon
+  real,dimension(nx_loc,ny_loc), intent(in) :: orog_lake_loc ! orographie a la surface des lacs
+  integer,dimension(80,2),intent(in) :: inclake_loc ! 80 voisins (+4 -4) en i,j
+  integer, dimension(nx_loc,ny_loc),intent(out) :: mask_extlake_loc ! extension du lac sous la glace
+  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)
+
+  integer, dimension(8) :: ipinctab,jpinctab ! indice des points a scanner autour du point ii jj
+  integer, dimension(nx_loc,ny_loc) :: clake ! point lac avec voisin point englace
+  integer :: i,j,k
 
 ! initialisation
-clake(:,:)=0
-orog_extlake_loc(:,:)=0.
-mask_extlake_loc(:,:)=9999
+  clake(:,:)=0
+  orog_extlake_loc(:,:)=0.
+  mask_extlake_loc(:,:)=9999
 
 ! recherche des points lacs en bordure de calotte => voisin pas lac englace et terre
-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. &
-  ((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. &
-  ((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. &
-  ((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))))
-  clake(:,:) = 1
-endwhere
+  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. &
+    ((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. &
+    ((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. &
+    ((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))))
+    clake(:,:) = 1
+  endwhere
 
 
 
 ! point terre et englace : mk=0 et H>0
-do j=1,ny
-  do i=1,nx
-    if (clake(i,j).EQ.1) then
-      do k=1,80
-        ipinctab(k) = max(1,min(nx,i + inclake_loc(k,1)))
-        jpinctab(k) = max(1,min(ny,j + inclake_loc(k,2)))
-        mask_extlake_loc(ipinctab(k),jpinctab(k)) = 1
-!        mask_extlake(ipinctab(k),jpinctab(k))=min(abs(inclake(k,1)),abs(inclake(k,2)),mask_extlake(ipinctab(k),jpinctab(k)))
-!        mask_extlake(ipinctab(k),jpinctab(k))=min(minval(abs(inclake(k,:),mask=inclake(k,:).GT.0)),mask_extlake(ipinctab(k),jpinctab(k)))
-        orog_extlake_loc(ipinctab(k),jpinctab(k))=orog_lake_loc(i,j)
-      enddo
-    endif
-  enddo
-enddo
-
-!~ print*,'!!!!!!! extlake !!!!!!'
-!~ print*, '50,138 lake,sealevel_2D,ice,mk,clake,mask_extlake,orog_extlake '
-!~ print*,lake_loc(50,138),ice_loc(50,138),mk_loc(50,138),mask_extlake_loc(50,138),orog_extlake_loc(50,138),clake(50,138)
-!~ print*,'65,166'
-!~ print*,lake_loc(65,166),ice_loc(65,166),mk_loc(65,166),mask_extlake_loc(65,166),orog_extlake_loc(65,166),clake(65,166)
-!~ print*,'64,166'
-!~ print*,lake_loc(64,166),ice_loc(64,166),mk_loc(64,166),mask_extlake_loc(64,166),orog_extlake_loc(64,166),clake(64,166)
-!~ print*,'66,166'
-!~ print*,lake_loc(66,166),ice_loc(66,166),mk_loc(66,166),mask_extlake_loc(66,166),orog_extlake_loc(66,166),clake(66,166)
-!~ print*,'65,165'
-!~ print*,lake_loc(65,165),ice_loc(65,165),mk_loc(65,165),mask_extlake_loc(65,165),orog_extlake_loc(65,165),clake(65,165)
-!~ print*,'65,167'
-!~ print*,lake_loc(65,167),ice_loc(65,167),mk_loc(65,167),mask_extlake_loc(65,167),orog_extlake_loc(65,167),clake(65,167)
+  do j=1,ny
+    do i=1,nx
+      if (clake(i,j).EQ.1) then
+        do k=1,80
+          ipinctab(k) = max(1,min(nx,i + inclake_loc(k,1)))
+          jpinctab(k) = max(1,min(ny,j + inclake_loc(k,2)))
+          mask_extlake_loc(ipinctab(k),jpinctab(k)) = 1
+          orog_extlake_loc(ipinctab(k),jpinctab(k))=orog_lake_loc(i,j)
+        enddo
+      endif
+    enddo
+  enddo
+
 end subroutine extlake