Changeset 110 for trunk


Ignore:
Timestamp:
05/24/17 17:46:02 (7 years ago)
Author:
aquiquet
Message:

Implementation of Schoof4Schoof parametrisation instead of Tsai4Schoof. Detection of floating points stuck inside grounded points (marais). Islands and marais are discarded for Schoof computation. flottab reverted to determin_tache and front (icebergs have to be put in calving for mass conservation).

Location:
trunk/SOURCES
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/3D-physique-gen_mod.f90

    r109 r110  
    489489  logical,dimension(nx,ny) :: new_flotmx !< pour signaler les points qui deviennent flottantmx entre 2 dtt  
    490490  logical,dimension(nx,ny) :: new_flotmy !< pour signaler les points qui deviennent flottantmy entre 2 dtt  
     491  logical,dimension(nx,ny) :: flot_marais  !< afq -- vrai si flottant et coince entre points poses 'o' 
    491492 
    492493  ! ===================== File id ========================================= 

  • trunk/SOURCES/Netcdf-routines/sortie_netcdf_GRISLI_mod.0.2-hassine.f90

    r107 r110  
    908908                   tab(:,:) = tot_water(:,:) 
    909909                endif 
    910                 if (itab.eq.59) then  
     910                if (itab.eq.59) then 
    911911                   tab(:,:) = gr_line_schoof(:,:) 
    912912                endif  
    … …  
    968968                            if (ilemy(i,j)) then 
    969969                               tab(i,j)=2 
     970                            else if (fleuvemy(i,j)) then 
     971                               tab(i,j)=5        ! actual grounded streams 
    970972                            else 
    971973                               tab(i,j)=1 

  • trunk/SOURCES/flottab2-0.7.f90

    r102 r110  
    5050 integer,dimension(0:n_ta_max)     :: nb_pts_tache !< indique le nombre de points par tache 
    5151 logical,dimension(0:n_ta_max)     :: iceberg      !< T si iceberg, F si calotte posee 
    52  
     52  
    5353 logical,dimension(nx,ny)        :: mask_tache_ij !< masque de travail  
    5454                                                  !< vrai pour toute la tache de i,j 
    … …  
    608608    !$OMP END PARALLEL 
    609609!print*,'H(90,179) flottab 3',H(90,179),ice(90,179),flot(90,179),Bsoc(90,179)+H(90,179)*ro/row -sealevel 
    610     !~ call DETERMIN_TACHE  
     610    call DETERMIN_TACHE  
    611611    !~  
    612612    !~ synchro :    if (isynchro.eq.1) then 
    … …  
    635635    !~ !!$ 
    636636    !~ !!$end if 
    637     !~  
    638     !~ !ice(:,:)=0  Attention, voir si ca marche toujours pour l'Antarctique et heminord ! 
    639     !~  
    640     !~ !On compte comme englacé uniquement les calottes dont une partie est posée       
    641     !~ !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij) 
    642     !~ !$OMP DO 
    643     !~ do j=3,ny-2 
    644     !~    do i=3,nx-2 
    645     !~ test1:  if (.not.iceberg(table_out(i,j))) then ! on est pas sur un iceberg  
    646     !~          if (nb_pts_tache(table_out(i,j)).ge.1) then 
    647     !~             ice(i,j)=1 
     637     
     638    !~ !ice(:,:)=0  ! Attention, voir si ca marche toujours pour l'Antarctique et heminord ! 
     639     
     640    !On compte comme englacé uniquement les calottes dont une partie est posée       
     641    !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij) 
     642    !$OMP DO 
     643    do j=3,ny-2 
     644       do i=3,nx-2 
     645    test1:  if (.not.iceberg(table_out(i,j))) then ! on est pas sur un iceberg  
     646             if (nb_pts_tache(table_out(i,j)).ge.1) then 
     647                 ice(i,j)=1 
    648648    !~             if (nb_pts_tache(table_out(i,j)).le.10) then  ! les petites iles sont en sia 
    649649    !~ !    write(6,*) 'petite ile ',i,j 
    … …  
    657657    !~  
    658658    !~  
    659     !~ ! ici on est sur une tache non iceberg >= 5 points                        
    660     !~ ! on teste si la tache n'est pas completement ice stream 
    661     !~  
    662     !~ test2:       if (icetrim(table_out(i,j))) then   ! on a une ile d'ice stream 
    663     !~  
    664     !~    mask_tache_ij(:,:)=.false. 
    665     !~    mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j))       ! pour toute la tache 
    666     !~  
    667     !~      smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) 
    668     !~      smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) 
    669     !~      smax_i=smax_coord(1) 
    670     !~      smax_j=smax_coord(2) 
    671     !~  
     659    ! ici on est sur une tache non iceberg >= 5 points                        
     660    ! on teste si la tache n'est pas completement ice stream 
     661     
     662    test2:       if (icetrim(table_out(i,j))) then   ! on a une ile d'ice stream 
     663     
     664       mask_tache_ij(:,:)=.false. 
     665       mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j))       ! pour toute la tache 
     666     
     667         smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:)) 
     668         smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:)) 
     669         smax_i=smax_coord(1) 
     670         smax_j=smax_coord(2) 
     671     
    672672    !~ !!$               smax_i=0 ; smax_j=0 ; smax_=sealevel 
    673673    !~ !!$               do ii=3,nx-2 
    … …  
    682682    !~ !!$                  end do 
    683683    !~ !!$               end do 
    684     !~  
    685     !~  
    686     !~                gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. 
    687     !~                gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. 
    688     !~                flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. 
    689     !~                flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 
    690     !~   
    691     !~                if (Smax_.le.sealevel) then 
    692     !~                   write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 
    693     !~                   write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j 
    694     !~                end if 
    695     !~  
    696     !~             endif test2  
    697     !~          end if                                             ! endif deplace 
    698     !~           
    699     !~       else ! on est sur un iceberg                          !   test1 
    700     !~          ice(i,j)=0 
    701     !~          h(i,j)=1. 
    702     !~          surnet=H(i,j)*(1.-ro/row) 
    703     !~          S(i,j)=surnet+sealevel 
    704     !~          B(i,j)=S(i,j)-H(i,j) 
    705     !~           
    706     !~       endif test1 
    707     !~  
    708     !~  
    709     !~    end do 
    710     !~ end do 
    711     !~ !$OMP END DO 
    712     !~ !$OMP END PARALLEL 
     684     
     685     
     686                   gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false. 
     687                   gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false. 
     688                   flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false. 
     689                   flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false. 
     690      
     691                   if (Smax_.le.sealevel) then 
     692                      write(num_tracebug,*)'Attention, une ile avec la surface sous l eau' 
     693                      write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j 
     694                   end if 
     695     
     696                endif test2  
     697             end if                                             ! endif deplace 
     698              
     699          else ! on est sur un iceberg                          !   test1 
     700             ice(i,j)=0 
     701             h(i,j)=0. !1. afq, we should put everything in calving! 
     702             surnet=H(i,j)*(1.-ro/row) 
     703             S(i,j)=surnet+sealevel 
     704             B(i,j)=S(i,j)-H(i,j) 
     705              
     706          endif test1 
     707     
     708     
     709       end do 
     710    end do 
     711     
     712    !$OMP END DO 
     713    !$OMP END PARALLEL 
     714 
    713715    !~  
    714716    !~ !---------------------------------------------- 
    … …  
    738740    !read(5,*) 
    739741 
    740     call determin_front_tof ! version simplifiee 
     742    call determin_front 
     743    !call determin_front_tof ! version simplifiee 
     744 
     745    call determin_marais 
    741746     
    742747  end subroutine flottab 
    … …  
    12921297end subroutine determin_front_tof 
    12931298 
     1299 
     1300!> SUBROUTINE: determin_marais 
     1301!! afq -- Routine pour l'identification des "marais" 
     1302!! Un marais est une tache "shelf" entouré de points grounded 
     1303!! Copie sauvage de determin_tache, adapte au probleme du marais 
     1304!> 
     1305subroutine determin_marais 
     1306 
     1307!$ USE OMP_LIB 
     1308 
     1309implicit none 
     1310integer :: indice 
     1311integer :: label         ! no des taches rencontrées dans le mask 
     1312integer :: label_max     ! no temporaire maxi de tache rencontrées 
     1313!      integer :: mask_nb = 4 
     1314integer,parameter :: mask_nb = 2   ! version ou on ne compte pas les diagonales 
     1315integer :: vartemp       ! variable temporaire pour reorganiser compt 
     1316!     integer,dimension(mask_nb) :: mask 
     1317integer,dimension(mask_nb) :: mask 
     1318integer pm ! loop integer 
     1319 
     1320integer,dimension(nx,ny)        :: table_out_marais      !< pour les numeros des taches 
     1321integer,dimension(0:n_ta_max)     :: compt_marais        !< contient les equivalence entre les taches 
     1322integer,dimension(0:n_ta_max)     :: nb_pts_marais !< indique le nombre de points par tache 
     1323logical,dimension(0:n_ta_max)     :: marais      !< T si iceberg, F si calotte posee 
     1324  
     1325! 1-initialisation 
     1326!----------------- 
     1327label_max=1      ! numero de la tache, la premiere tache est notée 1 
     1328label=1 
     1329do i=1,n_ta_max 
     1330   compt_marais(i)=i 
     1331enddo 
     1332!$OMP PARALLEL 
     1333!$OMP WORKSHARE 
     1334table_out_marais(:,:) = 0 
     1335marais(:)  = .true. 
     1336nb_pts_marais(:) = 0 
     1337!$OMP END WORKSHARE 
     1338!$OMP END PARALLEL 
     1339 
     1340! 2-reperage des taches 
     1341!---------------------- 
     1342 
     1343do j=2,ny-1 
     1344  do i=2,nx-1 
     1345 
     1346 
     1347 
     1348     IF ((ice(i,j).ge.1).and.flot(i,j)) THEN ! on est sur la glace qui flotte-------------------! 
     1349     !IF ((ice(i,j).ge.1).and.flot(i,j).and.(H(i,j).gt.1.)) THEN ! on est sur la glace qui flotte-------------------! 
     1350                 
     1351        !write (*,*) "un point qui nous interesse!",i,j 
     1352         
     1353        if (((ice(i-1,j).ge.1).and.flot(i-1,j)).or.((ice(i,j-1).ge.1).and.flot(i,j-1))) then   !masque de 2 cases adjacentes 
     1354        !if (((ice(i-1,j).ge.1).and.flot(i-1,j).and.(H(i-1,j).gt.1.)).or.((ice(i,j-1).ge.1).and.flot(i,j-1).and.(H(i,j-1).gt.1.))) then   !masque de 2 cases adjacentes 
     1355           !      un des voisins est deja en glace 
     1356           mask(1) = table_out_marais(i-1,j) 
     1357           mask(2) = table_out_marais(i,j-1)  
     1358           label = label_max 
     1359 
     1360           !on determine la valeur de la tache minimun (>0) presente ds le masque 
     1361           do indice=1,mask_nb 
     1362              if (mask(indice).gt.0) label=min(label,mask(indice)) 
     1363           enddo 
     1364 
     1365           !on fixe la valeur de la tache voisine minimun au point etudie (via label) 
     1366           table_out_marais(i,j)=label 
     1367            
     1368           !si un des voisins n'est pas glace alors la tache n'est pas un marais 
     1369           do pm=-1,1,2 
     1370              if ( (ice(i+pm,j).eq.0) .or. (ice(i,j+pm).eq.0) ) then 
     1371                 marais(label)=.false. 
     1372              endif 
     1373           enddo 
     1374 
     1375           ! si 2 taches differentes sont dans le masque, il faut les identifier dans compt 
     1376           ! i.e. les plus grands numeros correspondent au plus petit  
     1377           ! on lui affecte le numero de la tache fondamentale avec un signe - 
     1378           ! pour indiquer le changement 
     1379 
     1380           do indice=1,mask_nb 
     1381              if(mask(indice).gt.label) then 
     1382                 compt_marais(mask(indice))=-label 
     1383              endif 
     1384           enddo 
     1385            
     1386            
     1387        else !aucun des voisins est une tache 
     1388           table_out_marais(i,j)= label_max 
     1389           compt_marais(label_max)=label_max 
     1390            
     1391           !si un des voisins n'est pas glace alors la tache n'est pas un marais 
     1392           do pm=-1,1,2 
     1393              if ( (ice(i+pm,j).eq.0) .or. (ice(i,j+pm).eq.0) ) then 
     1394                 marais(label_max)=.false. 
     1395              endif 
     1396           enddo 
     1397 
     1398            
     1399           label_max  = label_max+1 
     1400           if (label_max.gt.n_ta_max) print*,'trop de taches=',label_max  
     1401        endif 
     1402 
     1403 
     1404     else           !on est pas sur une tache---------------------------------------------- 
     1405        table_out_marais(i,j)=0    ! Pas necessaire (reecrit 0 sur 0)              
     1406     endif          !--------------------------------------------------------------------- 
     1407 
     1408 
     1409  enddo 
     1410enddo 
     1411 
     1412 
     1413 
     1414! On reorganise compt en ecrivant le numero de la tache fondamentale 
     1415! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité) 
     1416! On indique aussi le nb de point que contient chaque taches (nb_pts_tache) 
     1417 
     1418do indice=1,label_max 
     1419   vartemp = compt_marais(indice) 
     1420   if (compt_marais(indice).lt.0) then  
     1421      compt_marais(indice)= compt_marais(-vartemp) 
     1422      if (.not.marais(indice)) marais(-vartemp)=.false. 
     1423   endif     
     1424enddo 
     1425 
     1426!$OMP PARALLEL 
     1427!$OMP DO REDUCTION(+:nb_pts_tache) 
     1428do j=1,ny 
     1429   do i=1,nx 
     1430      if (table_out_marais(i,j).ne.0) then 
     1431         table_out_marais(i,j)=compt_marais(table_out_marais(i,j))   
     1432         nb_pts_marais(compt_marais(table_out_marais(i,j)))= nb_pts_marais(compt_marais(table_out_marais(i,j)))+1    
     1433      endif 
     1434   enddo 
     1435enddo 
     1436!$OMP END DO 
     1437!$OMP END PARALLEL 
     1438 
     1439do j=1,ny 
     1440   do i=1,nx 
     1441      flot_marais(i,j) = marais(table_out_marais(i,j)) 
     1442   enddo 
     1443enddo 
     1444 
     1445end subroutine determin_marais 
     1446 
    12941447end module flottab_mod 
    12951448 

  • trunk/SOURCES/furst_schoof_mod.f90

    r109 r110  
    88implicit none 
    99 
    10 real :: frot_coef 
     10real :: frot_coef ! afq: fixed in Tsai but = beta in Schoof (for m_weert = 1) 
    1111real :: alpha_flot 
    1212integer :: gr_select 
    1313real,dimension(nx,ny) :: back_force_x 
    1414real,dimension(nx,ny) :: back_force_y 
     15real, parameter :: m_weert = 1.0 ! afq: is this defined elsewhere in the model??? 
    1516 
    1617contains 
    … …  
    3334 
    3435!  frot = 0.035 ! f in the solid friction law, 0.035 in Ritz et al. 2001         
    35    
     36 
    3637  ! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman) 
    3738  alpha_flot= ro/row 
    … …  
    4748  real :: Hgl 
    4849  !  integer :: igl 
    49  
    50   real :: back_force_coef=0.1 !!!!!! TO DO !!!!!!!! 
    5150 
    5251  real :: phi_im2, phi_im1, phi_i, phi_ip1, phi_ip2, phi_ip3 
    … …  
    8483  county(:,:)=0. 
    8584  gr_line_schoof(:,:)=0 
    86    
     85 
    8786  archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel 
    8887 
    … …  
    9796           ! selon x (indice i) 
    9897            
    99            if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0.) then 
     98           if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) then 
    10099           !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West) 
    101100 
    … …  
    106105              if (gr_select.eq.1) then ! flux de Tsai 
    107106                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_x(i-1,j),phi_prescr) 
     107              else if (gr_select.eq.2) then ! flux de Schoof 
     108                 ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point. 
     109                 if (xpos .lt. -dx/2.) then 
     110                    frot_coef = betamx(i,j) 
     111                 else 
     112                    frot_coef = betamx(i+1,j) 
     113                 endif 
     114                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,back_force_x(i-1,j),phi_prescr) 
    108115              else 
    109                  print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE' 
     116                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
    110117                 stop 
    111118              endif 
    … …  
    120127 
    121128 
    122  
     129                 if (.not.(ilemx(i,j)).and..not.(ilemx(i-1,j))) then 
     130                     
    123131                 ! extrapolation pour avoir uxbar(i-1,j) qui sera ensuite prescrit dans call rempli_L2 
    124132 
    … …  
    129137                 imx_diag (i-1,j)   = 1 
    130138                 gr_line_schoof(i-1,j) = 1 
    131                                                                   
     139 
    132140                 ! extrapolation pour avoir uxbar(i,j) 
    133141 
    … …  
    144152                 !                 print* 
    145153 
     154                 end if 
     155               
    146156              else            !  GL a Est du i staggered,                          o centre, x stag 
    147157 
    … …  
    150160                 !flux        !     in             in           out  G         out           in 
    151161 
    152  
     162                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then 
    153163 
    154164                 ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    … …  
    159169                 countx   (i,j)   = countx(i,j)+1 
    160170                 imx_diag (i,j) = 1 
    161                 gr_line_schoof(i,j) = 1 
    162                                                                   
     171                gr_line_schoof(i,j) = 1 
     172 
    163173                 ! extrapolation pour avoir uxbar(i+1,j) 
    164174 
    … …  
    169179                 imx_diag (i+1,j)   = 1 
    170180                 gr_line_schoof(i+1,j) = 1 
     181 
     182                 end if 
     183 
    171184              end if 
    172185                                          
    173            else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0.) then 
     186           else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) then 
    174187           !else if (flot(i+1,j).and.Uxbar(i+1,j).GT.0.) then    ! grounding line between i,i+1    CAS EST (ecoulement vers grille Est) 
    175188 
    … …  
    179192              if (gr_select.eq.1) then ! flux de Tsai 
    180193                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_x(i+1,j),phi_prescr) 
     194              else if (gr_select.eq.2) then ! flux de Schoof 
     195                 ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point. 
     196                 if (xpos .lt. -dx/2.) then 
     197                    frot_coef = betamx(i,j) 
     198                 else 
     199                    frot_coef = betamx(i+1,j) 
     200                 endif 
     201                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,back_force_x(i+1,j),phi_prescr) 
    181202              else 
    182                  print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE' 
     203                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
    183204                 stop 
    184205              endif 
    … …  
    190211                 !flux        !                    in           out        G  out            in 
    191212 
    192  
     213                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then 
     214                     
    193215                 ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    194216 
    … …  
    198220                 countx   (i,j)   = countx(i,j)+1 
    199221                 imx_diag (i,j) = 1 
    200                 gr_line_schoof(i,j) = 1 
    201                                                                   
     222                gr_line_schoof(i,j) = 1 
     223 
    202224                 ! extrapolation pour avoir uxbar(i+1,j) 
    203225 
    … …  
    207229                 countx   (i+1,j)   = countx(i+1,j)+1 
    208230                 imx_diag (i+1,j) = 1 
    209                  gr_line_schoof(i+1,j) = 1 
    210                                                                   
     231                 gr_line_schoof(i+1,j) = 1 
     232 
     233                 end if 
     234               
    211235              else   !  GL a west du i staggered,                          o centre, x stag 
    212236 
    … …  
    215239                 !flux        !                     in           out  G        out            in 
    216240 
    217  
     241                 if (.not.(ilemx(i+1,j)).and..not.(ilemx(i+2,j))) then 
     242                     
    218243                 ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    219244 
    … …  
    223248                 countx   (i+1,j)   = countx(i+1,j)+1 
    224249                 imx_diag (i+1,j) = 1 
    225                 gr_line_schoof(i+1,j) = 1 
     250                gr_line_schoof(i+1,j) = 1 
    226251 
    227252                 ! extrapolation pour avoir uxbar(i+1,j) 
    … …  
    233258                 imx_diag (i+2,j) = 1 
    234259                 gr_line_schoof(i+2,j) = 1 
     260 
     261                 end if 
     262               
    235263              end if 
    236264 
    … …  
    239267           !******************************************************************************* 
    240268           ! selon y (indice j) 
    241            if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0.) then  
     269           if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) then  
    242270           !if (flot(i,j-1).and.Uybar(i,j).LT.0.) then    ! grounding line between j-1,j    CAS SUD (ecoulement vers grille SUD) 
    243271 
    … …  
    248276              if (gr_select.eq.1) then ! flux de Tsai 
    249277                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_y(i,j-1),phi_prescr) 
     278              else if (gr_select.eq.2) then ! flux de Schoof 
     279                 ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point. 
     280                 if (ypos .lt. -dy/2.) then 
     281                    frot_coef = betamy(i,j) 
     282                 else 
     283                    frot_coef = betamy(i,j+1) 
     284                 endif 
     285                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,back_force_y(i,j-1),phi_prescr) 
    250286              else 
    251                  print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE' 
     287                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
    252288                 stop 
    253289              endif 
    … …  
    260296                 !flux        !     in            out         G out            in 
    261297 
     298                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j-1))) then 
    262299 
    263300                 ! extrapolation pour avoir uybar(i,j-1) qui sera ensuite prescrit dans call rempli_L2 
    … …  
    279316                 gr_line_schoof(i,j) = 1 
    280317 
     318 
     319                 end if 
     320               
    281321              else            !  GL au nord du j staggered,                          o centre, x stag 
    282322 
    … …  
    285325                 !flux        !     in             in           out  G         out           in 
    286326 
    287  
     327                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then 
     328                     
    288329                 ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    289330 
    … …  
    303344                 imy_diag (i,j+1) = 1 
    304345                 gr_line_schoof(i,j+1) = 1 
     346 
     347                 end if 
     348               
    305349              end if 
    306350               
    307            else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0.) then 
     351           else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) then 
    308352           !else if (flot(i,j+1).and.Uybar(i,j+1).GT.0.) then    ! grounding line between j,j+1    CAS NORD (ecoulement vers grille Nord) 
    309353 
    … …  
    313357              if (gr_select.eq.1) then ! flux de Tsai 
    314358                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_y(i,j+1),phi_prescr) 
     359              else if (gr_select.eq.2) then ! flux de Schoof 
     360                 ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point. 
     361                 if (ypos .lt. -dy/2.) then 
     362                    frot_coef = betamy(i,j) 
     363                 else 
     364                    frot_coef = betamy(i,j+1) 
     365                 endif 
     366                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,back_force_y(i,j+1),phi_prescr) 
    315367              else 
    316                  print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE' 
     368                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE' 
    317369                 stop 
    318370              endif 
    … …  
    324376                 !flux        !                    in           out        G  out            in 
    325377 
     378                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then 
    326379 
    327380                 ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    … …  
    343396                 gr_line_schoof(i,j+1) = 1 
    344397 
     398                 end if 
     399               
    345400              else   !  GL au nord du j staggered,                          o centre, x stag 
    346401 
    … …  
    349404                 !flux        !                     in           out  G        out            in 
    350405 
     406                 if (.not.(ilemy(i,j+1)).and..not.(ilemy(i,j+2))) then 
    351407 
    352408                 ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2 
    … …  
    357413                 county   (i,j+1)   = county(i,j+1)+1 
    358414                 imy_diag (i,j+1) = 1 
    359                 gr_line_schoof(i,j+1) = 1 
    360                                                                   
     415                gr_line_schoof(i,j+1) = 1 
     416 
    361417                 ! extrapolation pour avoir uybar(i,j+1) 
    362418 
    … …  
    367423                 imy_diag (i,j+2) = 1 
    368424                 gr_line_schoof(i,j+2) = 1 
    369                   
     425 
     426                 end if 
     427 
    370428              end if 
    371429 
    … …  
    376434 
    377435  ! afq -- if we discard the points with multiple fluxes coming, uncomment: 
    378   !where (countx(:,:).ge.2) 
    379   !   uxbartemp(:,:)=uxbar(i,j)*countx(:,:) 
    380   !   imx_diag(:,:) = imx_diag_in(:,:) 
    381   !end where 
    382   !where (county(:,:).ge.2) 
    383   !   uybartemp(:,:)=uybar(:,:)*county(:,:) 
    384   !   imy_diag(:,:) = imx_diag_in(:,:) 
    385   !end where 
     436  where (countx(:,:).ge.2) 
     437     uxbartemp(:,:)=uxbar(i,j)*countx(:,:) 
     438     imx_diag(:,:) = imx_diag_in(:,:) 
     439  end where 
     440  where (county(:,:).ge.2) 
     441     uybartemp(:,:)=uybar(:,:)*county(:,:) 
     442     imy_diag(:,:) = imx_diag_in(:,:) 
     443  end where 
    386444  ! afq -- 
    387445   
    … …  
    501559real,intent(out) :: phi_schoof 
    502560 
    503  
    504 phi_schoof = ((A_mean * rog**(n_glen+1) * (1 - alpha)**n_glen) / (4**n_glen *C_frot)) **(1./(1.+1./m_weert)) 
     561! afq: in standard GRISLI (19/05/17), tob= -beta ub => m_weert = 1 and C_frot = beta 
     562 
     563 
     564phi_schoof = (((A_mean * rog**(n_glen+1) * (1 - alpha)**n_glen) / (4**n_glen *C_frot)) **(1./(1.+1./m_weert)))*back_force_coef 
    505565phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert)) 
    506566 
Note: See TracChangeset for help on using the changeset viewer.