- Timestamp:
- 10/14/16 14:45:46 (8 years ago)
- Location:
- branches/iLoveclim/SOURCES
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/iLoveclim/SOURCES/3D-physique-gen_mod.f90
r89 r90 251 251 !hassine 252 252 real, dimension(nx,ny) :: beta_centre !< beta on major node (average) 253 real (kind = 8),dimension(nx,ny) :: calvingGRIS 253 254 real (kind = 8),dimension(nx,ny),target :: calvin_GRIS ! calving flux[m3/s] given to CLIO 254 255 real, dimension(nx,ny) :: coef_drag !< coefficient de la loi de friction non lineaire : depend de la valeur de alpha_drag -
branches/iLoveclim/SOURCES/Hemin40_files/output_hemin40_mod.f90
r77 r90 27 27 !real :: pdfmean ! moyenne PDF 28 28 logical,dimension(nx,ny,13) :: mask_cal ! masque regions calotte 29 REAL, dimension(nx,ny) :: old_H_dtt, old_H_new ! Epaisseur de glace au pas de temps precedent 30 real, dimension(nx,ny) :: runof_oc_sngl 29 31 integer, dimension(nx,ny) :: write_mask 30 32 CONTAINS … … 36 38 mask_cal(:,:,:)=.false. 37 39 mask_cal(:,:,1)=.true. 40 old_H_dtt(:,:)=H(:,:) ! afq missing in mab? 38 41 ! calcul d'un masque pour les regions des calottes 39 42 do j=1,ny … … 178 181 ! REAL ABLA(NX,NY) 179 182 REAL DELTAVOL 183 REAL, dimension(nx,ny,13) :: delta_H_dtt 184 REAL, dimension(nx,ny) :: delta_H_test 180 185 181 186 … … 207 212 TBM(kk) = 0. 208 213 ITJJA(kk) = 0. 214 delta_H_dtt(:,:,kk) = 0. 215 runof_oc(:,:) = 0. 209 216 210 217 ! nouveau tof mai 2009 … … 231 238 232 239 240 where (mask_cal(:,:,1).and.(H(:,:).gt.2.).and..not.flot(:,:)) 241 delta_H_dtt(:,:,1) = H(:,:) 242 elsewhere 243 delta_H_dtt(:,:,1) = 0d0 244 endwhere 245 old_H_new(:,:) = delta_H_dtt(:,:,1) 246 delta_H_dtt(:,:,1) = old_H_dtt(:,:) - delta_H_dtt(:,:,1) 247 248 runof_oc(:,:) = DBLE(delta_H_dtt(:,:,1))*DX*DY 249 250 runof_oc_sngl(:,:) = delta_H_dtt(:,:,1)*DX*DY 251 252 253 old_H_dtt(:,:) = old_H_new(:,:) 254 255 calvin_GRIS(:,:) = calvingGRIS(:,:) 256 calvingGRIS(:,:) =0. 257 233 258 !!$ 234 259 !!$ -
branches/iLoveclim/SOURCES/calving_frange.f90
r4 r90 24 24 real :: hcoup_max ! epaisseur max 25 25 real :: hcoup_min=100. ! epaisseur min 26 real :: mass_total 26 27 integer :: meth_hcoup ! pour avoir hcoup dépendant du coefbmshelf 27 28 integer :: ifrange ! 0 pas de traitement particulier pres du bord, 1 -> franges … … 31 32 logical :: avalw,avale,avals,avaln,interieur 32 33 34 integer :: testH ! somme des flux entrants en un point ij 35 33 36 contains 34 37 … … 40 43 41 44 45 mass_total = 0d0 46 calvingGRIS=0d0 42 47 calv(:,:)=0. 43 48 … … 72 77 !--------------------------------------------------------------------------------------- 73 78 subroutine calving 74 79 80 real :: calvtest 81 75 82 ! initialisation calving 76 83 calv(:,:)=0. 84 mass_total=0d0 85 calvtest=0. 86 mass_total=0d0 77 87 78 88 ! calcul du dhdt lagrangien … … 151 161 ! ifrange = 7 calving drastique sauf si coefbmshelf < 0.5 152 162 163 !cdc 164 ! ifrange=44 idem 4 mais avec test sur i et j couple // afq: probably useless 165 ! ifrange=8 C. Dumas : coupure si alimentation ne permet pas H > Hcoup 153 166 154 167 if (ifrange.eq.1) then ! Rajout vince : si un point voisin est pose -> test vrai. … … 315 328 .and.(coefbmshelf.lt.1.)) 316 329 317 330 else if (ifrange.eq.44) then ! idem 4 mais avec test sur i et j couple 331 332 bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.) 333 334 testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup 335 .and. (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) & 336 .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) ! 337 338 testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup 339 .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) & 340 .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) ! 341 342 testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup 343 .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))& 344 .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) ! 345 346 testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup 347 .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))& 348 .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) ! 349 350 else if (ifrange.eq.8) then ! test C Dumas : si somme flux entrant permet H > Hcoup : vrai 351 ! bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.) 352 353 testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)) & ! voisin (i-1,j) amont et > hcoup 354 ! .and. (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) & 355 .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) ! 356 357 testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.)) & ! voisin (i+1,j) amont et > hcoup 358 ! .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) & 359 .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) ! 360 361 testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)) & ! voisin (i,j-1) amont et > hcoup 362 ! .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))& 363 .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) ! 364 365 testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.)) & ! voisin (i,j+1) amont et > hcoup 366 ! .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))& 367 .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) ! 368 369 testH=0. 370 if (testmij) testH=testH + abs(uxbar(i,j))*H(i-1,j) 371 if (testpij) testH=testH + abs(uxbar(i+1,j))*H(i+1,j) 372 if (testimj) testH=testH + abs(uybar(i,j))*H(i,j-1) 373 if (testipj) testH=testH + abs(uybar(i,j+1))*H(i,j+1) 374 375 if (testH.gt.Hcoup) then 376 testmij=.true. 377 else 378 testmij=.false. 379 testpij=.false. 380 testimj=.false. 381 testipj=.false. 382 endif 383 318 384 endif 319 385 … … 326 392 iin3=max(1,i-3) 327 393 328 avalw=((hmhc(i-1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) & 329 .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i,j).le.0.) 330 394 !cdc avalw=((hmhc(i-1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) & 395 !cdc .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i,j).le.0.) 396 397 avalw=((hmhc(i-1,j).gt.0..and.flot(i-1,j)).or.(hmhc(iin2,j).gt.0..and.flot(iin2,j)) & 398 .or.(hmhc(iin3,j).gt.0..and.flot(iin3,j))).and.(uxbar(i,j).le.0.) 331 399 332 400 ! 1 des 3 voisins aval > hcoup du cote est (i+1) … … 334 402 iin3=min(i+3,nx) 335 403 336 avale=((hmhc(i+1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) & 337 .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i+1,j).ge.0) 338 404 !cdc avale=((hmhc(i+1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) & 405 !cdc .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i+1,j).ge.0) 406 407 avale=((hmhc(i+1,j).gt.0..and.flot(i+1,j)).or.(hmhc(iin2,j).gt.0..and.flot(iin2,j)) & 408 .or.(hmhc(iin3,j).gt.0..and.flot(iin3,j))).and.(uxbar(i+1,j).ge.0) 409 339 410 ! 1 des 3 voisins aval > hcoup du cote sud (j-1) 340 411 jin2=max(1,j-2) 341 412 jin3=max(1,j-3) 342 413 343 avals=((hmhc(i,j-1).gt.0.).or.(hmhc(i,jin2).gt.0.) &344 .or.(hmhc(i,jin3).gt.0.)).and.(uybar(i,j).le.0.)414 !cdc avals=((hmhc(i,j-1).gt.0.).or.(hmhc(i,jin2).gt.0.) & 415 !cdc .or.(hmhc(i,jin3).gt.0.)).and.(uybar(i,j).le.0.) 345 416 417 avals=((hmhc(i,j-1).gt.0..and.flot(i,j-1)).or.(hmhc(i,jin2).gt.0..and.flot(i,jin2)) & 418 .or.(hmhc(i,jin3).gt.0..and.flot(i,jin3))).and.(uybar(i,j).le.0.) 419 346 420 ! 1 des 3 voisins aval > hcoup du cote nord (j-1) 347 421 jin2=min(j+2,ny) 348 422 jin3=min(j+3,ny) 349 423 350 avaln=((hmhc(i,j+1).gt.0.).or.(hmhc(i,jin2).gt.0.) & 351 .or.(hmhc(i,jin2).gt.0.)) .and.(uybar(i,j+1).ge.0.) 424 !cdc avaln=((hmhc(i,j+1).gt.0.).or.(hmhc(i,jin2).gt.0.) & 425 !cdc .or.(hmhc(i,jin2).gt.0.)) .and.(uybar(i,j+1).ge.0.) 426 427 avaln=((hmhc(i,j+1).gt.0..and.flot(i,j+1)).or.(hmhc(i,jin2).gt.0..and.flot(i,jin2)) & 428 .or.(hmhc(i,jin2).gt.0..and.flot(i,jin2))) .and.(uybar(i,j+1).ge.0.) 352 429 353 430 interieur=(avalw.or.avale).and.(avals.or.avaln) … … 355 432 356 433 357 if ((.not.(testmij.or.testpij.or.testimj.or.testipj)) & ! pas suffisament alimente 358 .and.(.not.interieur)) then ! et pas interieur 434 !cdc if ((.not.(testmij.or.testpij.or.testimj.or.testipj)) & ! pas suffisament alimente 435 !cdc .and.(.not.interieur)) then ! et pas interieur 436 437 !cdc pour ne pas faire de calving si H > 350m : if ((.not.(testmij.or.testpij.or.testimj.or.testipj)).and.(H(i,j).LT.350)) then 438 ! mab: if none of the neighbours doesn't provide ice for the considered point, 439 ! mab: then the ice is cut at this point and given to calv(i,j) 440 if ((.not.(testmij.or.testpij.or.testimj.or.testipj))) then 441 ! mab: new version so that grid points with h=1m (all points were calving takes place are afterwards 442 ! mab: set to a height of 1m) are not taken into account, provides a problem as soon as you start 443 ! mab: coupling GRISLI to the iceberg modell 444 if(h(i,j).gt.1) then 359 445 calv(i,j)=-h(i,j) 360 H(i,j)=1 446 ! mab: copy the calving in the right format for coupling 447 ! calvin_GRIS = volume flux per year (m3/a) 448 calvin_GRIS(i,j)=calvin_GRIS(i,j)+(-calv(i,j)*dx*dy) 449 ! calvin_GRIS(i,j)=-calv(i,j)*40000.*40000. 450 ! write(15,*) 'calvheight:',h(i,j),calv(i,j),i,j,dtt,dt,calvin_GRIS(i,j) 361 451 endif 452 H(i,j)=1 453 endif 362 454 363 455 end if ifint … … 365 457 end do 366 458 end do 459 460 calvtest=sum(-calv(:,:)*dx*dy) 461 mass_total=sum(calvingGRIS(:,:)) 462 463 calvtest = 0. 464 465 367 466 end subroutine calving 368 467 !------------------------------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.