Changeset 96 for trunk/SOURCES


Ignore:
Timestamp:
12/09/16 10:58:15 (8 years ago)
Author:
dumas
Message:

Schoof update : back force is now calculated by 2 calls in diagnoshelf

Location:
trunk/SOURCES
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/New-remplimat/diagno-L2_mod.f90

    r93 r96  
    1212real, dimension(nx,ny) :: uxb1 
    1313real, dimension(nx,ny) :: uyb1 
     14 
     15real, dimension(nx,ny) :: uxb1ramollo 
     16real, dimension(nx,ny) :: uyb1ramollo 
     17real, dimension(nx,ny) :: pvi_keep 
     18 
    1419 
    1520!cdc transfere dans module3d pour compatibilite avec furst_schoof_mod 
    … …  
    178183  call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), &  
    179184       uxb1(nxd1:nxd2,nyd1:nyd2),uyb1(nxd1:nxd2,nyd1:nyd2),  & 
    180        imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno)     
     185       imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno) 
     186        
     187  if (Schoof.eq.1) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo 
     188    pvi_keep(:,:)=pvi(:,:)        
     189                where (flot(:,:)) 
     190                                pvi(:,:)=pvimin 
     191                endwhere                         
     192             
     193                call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), &  
     194       uxb1ramollo(nxd1:nxd2,nyd1:nyd2),uyb1ramollo(nxd1:nxd2,nyd1:nyd2),  & 
     195       imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno) 
     196    pvi(:,:)=pvi_keep(:,:)  
     197    where ((uxb1ramollo(:,:)**2 + uyb1ramollo(:,:)**2).GT.0.) 
     198                        back_force(:,:) = sqrt(uxb1(:,:)**2+ uyb1(:,:)**2) / sqrt(uxb1ramollo(:,:)**2 + uyb1ramollo(:,:)**2) 
     199                elsewhere 
     200                        back_force(:,:)=1. 
     201                endwhere 
     202  endif 
     203   
     204!~   do j=1,ny 
     205!~              do i=1,nx 
     206!~                if (sqrt(uxb1(i,j)**2+ uyb1(i,j)**2).gt.0..and..not.flot(i,j)) then 
     207!~                              write(1034,*) sqrt(uxb1(i,j)**2+ uyb1(i,j)**2) / sqrt(uxb1ramollo(i,j)**2 + uyb1ramollo(i,j)**2) 
     208!~                else 
     209!~                              write(1034,*) 1. 
     210!~                      endif 
     211!~              enddo            
     212!~      enddo 
     213         
     214!~   print*,'apres calcul rempli_L2' 
     215!~   read(*,*)       
    181216 
    182217  ! Dans rempli_L2 

  • trunk/SOURCES/furst_schoof_mod.f90

    r94 r96  
    1111real :: alpha_flot 
    1212integer :: gr_select 
     13real,dimension(nx,ny) :: back_force 
    1314 
    1415contains 
    … …  
    3435  ! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman) 
    3536  alpha_flot= ro/row 
     37  back_force(:,:)=0.1 
    3638   
    3739end subroutine init_furst_schoof   
    … …  
    6466 
    6567  real,dimension(nx,ny) :: uxbartemp, uybartemp 
     68   
     69  real,dimension(nx,ny) :: archimtab 
    6670 
    6771  phi_prescr_tabx=0. 
    … …  
    7579  countx(:,:)=0. 
    7680  county(:,:)=0. 
     81 
     82        archimtab(:,:) = Bsoc(:,:)+H(:,:)*ro/row -sealevel 
    7783 
    7884  ! detection des noeuds grounding line 
    … …  
    8389           !           print*,'schoof gr line ij',i,j 
    8490           ! selon x (indice i) 
    85            if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West) 
     91            
     92           if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0.) then 
     93           !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West) 
    8694 
    8795              call  calc_xgl4schoof(-dx,alpha_flot,H(i,j),Bsoc(i,j),H(i-1,j),Bsoc(i-1,j),Sealevel,xpos,Hgl) 
    … …  
    9098 
    9199              if (gr_select.eq.1) then ! flux de Tsai 
    92                  call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_coef,phi_prescr) 
     100                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force(i-1,j),phi_prescr) 
    93101              else 
    94102                 print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE' 
    … …  
    151159                 imx_diag (i+1,j)   = 1 
    152160              end if 
    153  
    154            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) 
     161                                          
     162                                         else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0.) then 
     163           !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) 
    155164 
    156165              call  calc_xgl4schoof(dx,alpha_flot,H(i,j),Bsoc(i,j),H(i+1,j),Bsoc(i+1,j),Sealevel,xpos,Hgl) 
    … …  
    158167              Hglx_tab(i,j)=Hgl 
    159168              if (gr_select.eq.1) then ! flux de Tsai 
    160                  call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_coef,phi_prescr) 
     169                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force(i+1,j),phi_prescr) 
    161170              else 
    162171                 print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE' 
    … …  
    215224           !******************************************************************************* 
    216225           ! selon y (indice j) 
    217            if (flot(i,j-1).and.Uybar(i,j).LT.0.) then    ! grounding line between j-1,j    CAS SUD (ecoulement vers grille SUD) 
     226           if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0.) then  
     227           !if (flot(i,j-1).and.Uybar(i,j).LT.0.) then    ! grounding line between j-1,j    CAS SUD (ecoulement vers grille SUD) 
    218228 
    219229 
    … …  
    222232              Hgly_tab(i,j)=Hgl 
    223233              if (gr_select.eq.1) then ! flux de Tsai 
    224                  call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_coef,phi_prescr) 
     234                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force(i,j-1),phi_prescr) 
    225235              else 
    226236                 print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE' 
    … …  
    275285                 imy_diag (i,j+1) = 1 
    276286              end if 
    277  
    278            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) 
     287               
     288           else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0.) then 
     289           !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) 
    279290 
    280291              call  calc_xgl4schoof(dy,alpha_flot,H(i,j),Bsoc(i,j),H(i,j+1),Bsoc(i,j+1),Sealevel,ypos,Hgl) 
    … …  
    282293              Hgly_tab(i,j)=Hgl 
    283294              if (gr_select.eq.1) then ! flux de Tsai 
    284                  call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_coef,phi_prescr) 
     295                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force(i,j+1),phi_prescr) 
    285296              else 
    286297                 print*,'ATTENTION FLUX AUTRE QUE TSAI NON IMPLEMENTE' 
    … …  
    392403  real                   ::    Cgl       ! variable de travail   
    393404 
    394   Cgl   = (-alpha * (H_1-H_0) - (B_1-B_0)) ! -(rho/rhow (H1-H0) + (B1-B0)) 
     405!~   Cgl   = (-alpha * (H_1-H_0) - (B_1-B_0)) ! -(rho/rhow (H1-H0) + (B1-B0)) 
     406 
     407!~   if (abs(Cgl).gt.1.e-5) then    
     408!~      xpos    = (B_0 + alpha * H_0 - SL) / Cgl  
     409!~   else 
     410!~      xpos    = 1. - 1./abs(dy)   ! verifier 
     411!~   end if 
     412   
     413!~   if (xpos.LT.0..OR.xpos.GT.1.) then 
     414!~              print*,'calc_xgl4schoof : xpos value error i,j=',i,j 
     415!~              print*, 'xpos,Cgl', xpos,Cgl 
     416!~              print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL 
     417!~              stop 
     418!~      endif 
     419                 
     420  Cgl   = (B_1-B_0) + alpha * (H_1-H_0) 
    395421 
    396422  if (abs(Cgl).gt.1.e-5) then    
    397      xpos    = (B_0 + alpha * H_0 - SL) / Cgl  
     423     xpos    = (SL - (B_0 + alpha * H_0)) / Cgl  
    398424  else 
    399      xpos    = 1. - 1./abs(dy)   ! verifier 
     425     xpos    = 0.5  ! verifier 
    400426  end if 
    401427   
     428!  print* 
     429!  print*, 'xpos', xpos 
     430!       print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL 
     431!       print*,'B_1, H_1', B_1, H_1 
     432 
     433         
    402434  if (xpos.LT.0..OR.xpos.GT.1.) then 
    403435                print*,'calc_xgl4schoof : xpos value error i,j=',i,j 
     436                print*, 'xpos,Cgl', xpos,Cgl 
     437                print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL 
     438                print*,'archim=',B_1+H_1*alpha - SL 
    404439                stop 
    405         endif    
    406  
     440        endif 
     441         
    407442  xpos = xpos * dy 
    408443   
Note: See TracChangeset for help on using the changeset viewer.