Changeset 375 for branches/GRISLIv3


Ignore:
Timestamp:
03/20/23 15:55:43 (16 months ago)
Author:
dumas
Message:

use only in subroutine SIA_velocities

Location:
branches/GRISLIv3/SOURCES
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/3D-physique-gen_mod.f90

    r333 r375  
    350350  real,dimension(nx,ny) :: W0           !< enfoncement du socle a l'equilibre isostatique 
    351351  real,dimension(nx,ny) :: W1           !< enfoncement du socle courant 
    352   real,dimension(nx,ny) :: XX           !< work array 
    353352  real,dimension(nx,ny) :: XLONG        !< longitude 
    354353  real,dimension(nx,ny) :: YLAT         !< latitude 

  • branches/GRISLIv3/SOURCES/initial2-0.4.f90

    r237 r375  
    5252  UYBAR(:,:)=0. 
    5353  VBAR(:,:)=0. 
    54   XX(:,:)=0.  
    5554  IBASE(:,:)=1 
    5655  TPMP(:,:,1)=0 

  • branches/GRISLIv3/SOURCES/velocities-polyn-0.3.f90

    r372 r375  
    1515! ========================================================================== 
    1616! 
    17       subroutine SIA_velocities() 
     17subroutine SIA_velocities 
    1818! 
    1919! 
    … …  
    2323!     Attention la vitesse verticale est "super reduite" 
    2424! ========================================================================== 
    25 !$ USE OMP_LIB 
    26 use module3d_phy 
    27 !use deform_declar  
    28 USE module_choix 
     25  !$ USE OMP_LIB 
     26  use runparam, only : itracebug       
     27  use module3d_phy, only: nx,ny,nz,dx,flotmx,flotmy,ux,uy,sux,suy,sdx,sdy,ubx,uby,& 
     28       hmx,hmy,uxbar,uybar,iglen,cde,divu,uzr,bmelt,bm,flot,front,flgzmx,flgzmy,& 
     29       uzsdot,dtt,uzk,bdot,hdot 
     30  use deformation_mod_2lois, only: n1poly,n2poly,ddx,ddy,sa_mx,sa_my,s2a_mx,s2a_my 
    2931 
    30 implicit none 
     32  implicit none 
    3133 
    32 real, dimension(nx,ny) :: hdd 
     34  integer :: i,j,k 
     35  real, dimension(nx,ny) :: hdd 
     36  real, dimension(nx,ny) :: xx 
    3337 
    34 if (itracebug.eq.1)  call tracebug(' Entree dans routine SIA_velocity') 
     38  if (itracebug.eq.1)  call tracebug(' Entree dans routine SIA_velocity') 
    3539 
    36 !$OMP PARALLEL 
    37 !$OMP DO COLLAPSE(2) 
    38 do k=1,nz 
    39    do j=2,ny-1 
    40       do i=2,nx-1 
     40  !$OMP PARALLEL 
     41  xx(:,:)=0. 
     42  !$OMP DO COLLAPSE(2) 
     43  do k=1,nz 
     44     do j=2,ny-1 
     45        do i=2,nx-1 
    4146 
    42 !      flgzmx pour tous les points qui sont traites par l'equation elliptique 
    43 !      diagnoshelf 
     47           !      flgzmx pour tous les points qui sont traites par l'equation elliptique 
     48           !      diagnoshelf 
    4449 
    45 !     Selon X 
    46 !      ------------------------------------- 
    47 !         if (.not.flgzmx(i,j)) then      ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) 
    48 !            ux(i,j,k)=ddbx(i,j) 
    49 !            sux(i,j,k)=ddbx(i,j)*cde(k) 
     50           !     Selon X 
     51           !      ------------------------------------- 
     52           !         if (.not.flgzmx(i,j)) then      ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) 
     53           !            ux(i,j,k)=ddbx(i,j) 
     54           !            sux(i,j,k)=ddbx(i,j)*cde(k) 
    5055 
    51         if (.not.flotmx(i,j)) then        ! version aout 2010 (Cat)  
    52                                           ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) 
    53                                           ! et agit sur tous les points poses 
     56           if (.not.flotmx(i,j)) then        ! version aout 2010 (Cat)  
     57              ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) 
     58              ! et agit sur tous les points poses 
    5459 
    55             ux(i,j,k) = 0. 
    56             sux(i,j,k)= 0. 
     60              ux(i,j,k) = 0. 
     61              sux(i,j,k)= 0. 
    5762 
    58             do iglen=n1poly,n2poly 
    59                ux(i,j,k)  = ux(i,j,k)+ddx(i,j,iglen)*sa_mx(i,j,k,iglen) 
    60                sux(i,j,k) = sux(i,j,k)+ddx(i,j,iglen)*s2a_mx(i,j,k,iglen) 
    61             end do 
    62             ux(i,j,k)  = ux(i,j,k)*(-sdx(i,j)) 
    63             ux(i,j,k)  = ux(i,j,k) + ubx(i,j)    ! ajoute le glissement 
     63              do iglen=n1poly,n2poly 
     64                 ux(i,j,k)  = ux(i,j,k)+ddx(i,j,iglen)*sa_mx(i,j,k,iglen) 
     65                 sux(i,j,k) = sux(i,j,k)+ddx(i,j,iglen)*s2a_mx(i,j,k,iglen) 
     66              end do 
     67              ux(i,j,k)  = ux(i,j,k)*(-sdx(i,j)) 
     68              ux(i,j,k)  = ux(i,j,k) + ubx(i,j)    ! ajoute le glissement 
    6469 
    65             sux(i,j,k) = sux(i,j,k)*sdx(i,j)*hmx(i,j) 
    66             sux(i,j,k) = sux(i,j,k) - ubx(i,j) * cde(k) 
     70              sux(i,j,k) = sux(i,j,k)*sdx(i,j)*hmx(i,j) 
     71              sux(i,j,k) = sux(i,j,k) - ubx(i,j) * cde(k) 
    6772 
    68          else 
    69           sux(i,j,k) = (nz-k)*1./(nz-1.)*uxbar(i,j)*hmx(i,j) 
    70           ux(i,j,k)  = uxbar(i,j) 
    71          end if 
     73           else 
     74              sux(i,j,k) = (nz-k)*1./(nz-1.)*uxbar(i,j)*hmx(i,j) 
     75              ux(i,j,k)  = uxbar(i,j) 
     76           end if 
    7277 
    73 !     Selon Y 
    74 !      ------------------------------------- 
    75 !         if (.not.flgzmy(i,j)) then         
    76 !            uy(i,j,k)=ddby(i,j) 
    77 !            suy(i,j,k)=ddby(i,j)*cde(k) 
     78           !     Selon Y 
     79           !      ------------------------------------- 
     80           !         if (.not.flgzmy(i,j)) then         
     81           !            uy(i,j,k)=ddby(i,j) 
     82           !            suy(i,j,k)=ddby(i,j)*cde(k) 
    7883 
    79         if (.not.flotmy(i,j)) then        ! version aout 2010 (Cat)  
    80                                           ! utilise uby au lieu de ddby (compatibilite avec mix et spinup) 
    81                                           ! et agit sur tous les points poses 
     84           if (.not.flotmy(i,j)) then        ! version aout 2010 (Cat)  
     85              ! utilise uby au lieu de ddby (compatibilite avec mix et spinup) 
     86              ! et agit sur tous les points poses 
    8287 
    83             uy(i,j,k) = 0. 
    84             suy(i,j,k)= 0. 
     88              uy(i,j,k) = 0. 
     89              suy(i,j,k)= 0. 
    8590 
    8691 
    87             do iglen=n1poly,n2poly 
    88                uy(i,j,k)  = uy(i,j,k)+ddy(i,j,iglen)*sa_my(i,j,k,iglen) 
    89                suy(i,j,k) = suy(i,j,k)+ddy(i,j,iglen)*s2a_my(i,j,k,iglen) 
    90             end do 
    91             uy(i,j,k)  = uy(i,j,k)*(-sdy(i,j)) 
    92             uy(i,j,k)  = uy(i,j,k) + uby(i,j)    ! ajoute le glissement 
     92              do iglen=n1poly,n2poly 
     93                 uy(i,j,k)  = uy(i,j,k)+ddy(i,j,iglen)*sa_my(i,j,k,iglen) 
     94                 suy(i,j,k) = suy(i,j,k)+ddy(i,j,iglen)*s2a_my(i,j,k,iglen) 
     95              end do 
     96              uy(i,j,k)  = uy(i,j,k)*(-sdy(i,j)) 
     97              uy(i,j,k)  = uy(i,j,k) + uby(i,j)    ! ajoute le glissement 
    9398 
    94             suy(i,j,k) = suy(i,j,k)*sdy(i,j)*hmy(i,j) 
    95             suy(i,j,k) = suy(i,j,k) - uby(i,j) * cde(k) 
     99              suy(i,j,k) = suy(i,j,k)*sdy(i,j)*hmy(i,j) 
     100              suy(i,j,k) = suy(i,j,k) - uby(i,j) * cde(k) 
    96101 
    97          else 
    98           suy(i,j,k) = (nz-k)*1./(nz-1)*uybar(i,j)*hmy(i,j) 
    99           uy(i,j,k)=uybar(i,j) 
    100          endif 
     102           else 
     103              suy(i,j,k) = (nz-k)*1./(nz-1)*uybar(i,j)*hmy(i,j) 
     104              uy(i,j,k)=uybar(i,j) 
     105           endif 
    101106 
    102       end do 
    103    end do 
    104 end do 
    105 !$OMP END DO 
     107        end do 
     108     end do 
     109  end do 
     110  !$OMP END DO 
    106111 
    107 !       *************************** Z VELOCITIES ****************** 
     112  !       *************************** Z VELOCITIES ****************** 
    108113 
    109 !$OMP DO 
    110 do j=2,ny-1 
    111    do i=2,nx-1 
     114  !$OMP DO 
     115  do j=2,ny-1 
     116     do i=2,nx-1 
    112117 
    113       divu(i,j)=((uxbar(i+1,j)*hmx(i+1,j)-uxbar(i,j)*hmx(i,j))  &    ! formulation explicite 
    114            + (uybar(i,j+1)*hmy(i,j+1)-uybar(i,j)*hmy(i,j)))/dx 
     118        divu(i,j)=((uxbar(i+1,j)*hmx(i+1,j)-uxbar(i,j)*hmx(i,j))  &    ! formulation explicite 
     119             + (uybar(i,j+1)*hmy(i,j+1)-uybar(i,j)*hmy(i,j)))/dx 
    115120 
    116       uzr(i,j,nz)=bmelt(i,j) 
    117       xx(i,j)=uzr(i,j,1) 
    118       hdd(i,j)=bm(i,j)-bmelt(i,j)-divu(i,j)                          ! formulation explicite de dH/dt 
     121        uzr(i,j,nz)=bmelt(i,j) 
     122        xx(i,j)=uzr(i,j,1) 
     123        hdd(i,j)=bm(i,j)-bmelt(i,j)-divu(i,j)                          ! formulation explicite de dH/dt 
    119124 
    120 !     ------------------------------------------------------------------- 
    121 !      attention uzr contient maintenant : 
    122 !      uzr=uz +u ((1-e)*dH/dx+dB/dx)+ ((1-e)*dH/dt+dB/dt) 
     125        !     ------------------------------------------------------------------- 
     126        !      attention uzr contient maintenant : 
     127        !      uzr=uz +u ((1-e)*dH/dx+dB/dx)+ ((1-e)*dH/dt+dB/dt) 
    123128 
    124          do k=1,nz-1 
    125             if (.not.flot(i,j).and.(front(i,j).eq.4.).and..not.flgzmx(i,j) & 
    126                  .and..not.flgzmx(i+1,j).and..not.flgzmy(i,j) & 
    127                  .and..not.flgzmy(i,j+1)) then 
    128                uzr(i,j,k)=uzr(i,j,nz)-((sux(i+1,j,k)-sux(i,j,k)) & 
    129                 +(suy(i,j+1,k)-suy(i,j,k)))/dx 
    130                uzr(i,j,k)=uzr(i,j,k)+hdd(i,j)*cde(k) 
     129        do k=1,nz-1 
     130           if (.not.flot(i,j).and.(front(i,j).eq.4.).and..not.flgzmx(i,j) & 
     131                .and..not.flgzmx(i+1,j).and..not.flgzmy(i,j) & 
     132                .and..not.flgzmy(i,j+1)) then 
     133              uzr(i,j,k)=uzr(i,j,nz)-((sux(i+1,j,k)-sux(i,j,k)) & 
     134                   +(suy(i,j+1,k)-suy(i,j,k)))/dx 
     135              uzr(i,j,k)=uzr(i,j,k)+hdd(i,j)*cde(k) 
    131136 
    132             else 
     137           else 
    133138 
    134 !     dans l'ice-shelf la vitesse uzr est lineaire ... 
    135 !    ------------------------------------------------- 
     139              !     dans l'ice-shelf la vitesse uzr est lineaire ... 
     140              !    ------------------------------------------------- 
    136141 
    137 !     attention, dans les zones stream, le vitesse verticale n'est pas 
    138 !     lineaire si elle est voisine d'un point .not.flgzmx 
    139 !     le calcul se fait alors par differences finies 
     142              !     attention, dans les zones stream, le vitesse verticale n'est pas 
     143              !     lineaire si elle est voisine d'un point .not.flgzmx 
     144              !     le calcul se fait alors par differences finies 
    140145 
    141                uzr(i,j,k)= bm(i,j)+(k-1.)/(nz-1.)*(bmelt(i,j)-bm(i,j)) 
     146              uzr(i,j,k)= bm(i,j)+(k-1.)/(nz-1.)*(bmelt(i,j)-bm(i,j)) 
    142147           end if 
    143148 
    144149        end do 
    145          
     150 
    146151     end do 
    147  end do 
    148 !$OMP END DO 
     152  end do 
     153  !$OMP END DO 
    149154 
    150 !$OMP DO 
     155  !$OMP DO 
    151156  do j=2,ny-1 
    152157     do i=2,nx-1   
    … …  
    155160     end do 
    156161  end do 
    157 !$OMP END DO 
    158 !$OMP END PARALLEL 
     162  !$OMP END DO 
     163  !$OMP END PARALLEL 
    159164 
    160165end subroutine SIA_velocities 
Note: See TracChangeset for help on using the changeset viewer.