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