[4] | 1 | !> \file velocities-polyn-0.3.f90 |
---|
| 2 | !! Calcule des vitesses. |
---|
| 3 | !< |
---|
| 4 | |
---|
| 5 | !> SUBROUTINE: VELOCITIES() |
---|
| 6 | !! \author Catherine RITZ |
---|
| 7 | !! \date mars 2007 |
---|
| 8 | !! @note Cette routine permet de calculer les vitesses. |
---|
| 9 | !! @note Used modules: |
---|
| 10 | !! @note - use module3D_phy |
---|
| 11 | !! @note - use module_choix |
---|
| 12 | !! |
---|
| 13 | !< |
---|
| 14 | |
---|
| 15 | ! ========================================================================== |
---|
| 16 | ! |
---|
| 17 | subroutine SIA_velocities() |
---|
| 18 | ! |
---|
| 19 | ! |
---|
| 20 | ! Cat passage en f90 en mars 2007 |
---|
| 21 | ! moyenne sa_mx et sa_my |
---|
| 22 | ! |
---|
| 23 | ! Attention la vitesse verticale est "super reduite" |
---|
| 24 | ! ========================================================================== |
---|
[72] | 25 | !$ USE OMP_LIB |
---|
[4] | 26 | use module3d_phy |
---|
| 27 | !use deform_declar |
---|
| 28 | USE module_choix |
---|
| 29 | implicit none |
---|
| 30 | |
---|
| 31 | REAL hdd |
---|
| 32 | if (itracebug.eq.1) call tracebug(' Entree dans routine SIA_velocity') |
---|
| 33 | |
---|
[72] | 34 | !$OMP PARALLEL |
---|
| 35 | !$OMP DO COLLAPSE(2) |
---|
[4] | 36 | do k=1,nz |
---|
| 37 | do j=2,ny-1 |
---|
| 38 | do i=2,nx-1 |
---|
| 39 | |
---|
| 40 | ! flgzmx pour tous les points qui sont traites par l'equation elliptique |
---|
| 41 | ! diagnoshelf |
---|
| 42 | |
---|
| 43 | ! Selon X |
---|
| 44 | ! ------------------------------------- |
---|
| 45 | ! if (.not.flgzmx(i,j)) then ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) |
---|
| 46 | ! ux(i,j,k)=ddbx(i,j) |
---|
| 47 | ! sux(i,j,k)=ddbx(i,j)*cde(k) |
---|
| 48 | |
---|
| 49 | if (.not.flotmx(i,j)) then ! version aout 2010 (Cat) |
---|
| 50 | ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) |
---|
| 51 | ! et agit sur tous les points poses |
---|
| 52 | |
---|
| 53 | ux(i,j,k) = 0. |
---|
| 54 | sux(i,j,k)= 0. |
---|
| 55 | |
---|
| 56 | do iglen=n1poly,n2poly |
---|
| 57 | |
---|
| 58 | ux(i,j,k) = ux(i,j,k)+ddx(i,j,iglen)*sa_mx(i,j,k,iglen) |
---|
| 59 | sux(i,j,k) = sux(i,j,k)+ddx(i,j,iglen)*s2a_mx(i,j,k,iglen) |
---|
| 60 | end do |
---|
| 61 | ux(i,j,k) = ux(i,j,k)*(-sdx(i,j)) |
---|
| 62 | ux(i,j,k) = ux(i,j,k) + ubx(i,j) ! ajoute le glissement |
---|
| 63 | |
---|
| 64 | sux(i,j,k) = sux(i,j,k)*sdx(i,j)*hmx(i,j) |
---|
| 65 | sux(i,j,k) = sux(i,j,k) - ubx(i,j) * cde(k) |
---|
| 66 | |
---|
| 67 | else |
---|
| 68 | sux(i,j,k) = (nz-k)*1./(nz-1.)*uxbar(i,j)*hmx(i,j) |
---|
| 69 | ux(i,j,k) = uxbar(i,j) |
---|
| 70 | end if |
---|
| 71 | |
---|
| 72 | ! Selon Y |
---|
| 73 | ! ------------------------------------- |
---|
| 74 | ! if (.not.flgzmy(i,j)) then |
---|
| 75 | ! uy(i,j,k)=ddby(i,j) |
---|
| 76 | ! suy(i,j,k)=ddby(i,j)*cde(k) |
---|
| 77 | |
---|
| 78 | if (.not.flotmy(i,j)) then ! version aout 2010 (Cat) |
---|
| 79 | ! utilise uby au lieu de ddby (compatibilite avec mix et spinup) |
---|
| 80 | ! et agit sur tous les points poses |
---|
| 81 | |
---|
| 82 | uy(i,j,k) = 0. |
---|
| 83 | suy(i,j,k)= 0. |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | do iglen=n1poly,n2poly |
---|
| 87 | uy(i,j,k)=uy(i,j,k)+ddy(i,j,iglen)*sa_my(i,j,k,iglen) |
---|
| 88 | suy(i,j,k)=suy(i,j,k)+ddy(i,j,iglen)*s2a_my(i,j,k,iglen) |
---|
| 89 | end do |
---|
| 90 | uy(i,j,k) = uy(i,j,k)*(-sdy(i,j)) |
---|
| 91 | uy(i,j,k) = uy(i,j,k) + uby(i,j) ! ajoute le glissement |
---|
| 92 | |
---|
| 93 | suy(i,j,k) = suy(i,j,k)*sdy(i,j)*hmy(i,j) |
---|
| 94 | suy(i,j,k) = suy(i,j,k) - uby(i,j) * cde(k) |
---|
| 95 | |
---|
| 96 | else |
---|
| 97 | suy(i,j,k) = (nz-k)*1./(nz-1)*uybar(i,j)*hmy(i,j) |
---|
| 98 | uy(i,j,k)=uybar(i,j) |
---|
| 99 | endif |
---|
| 100 | |
---|
| 101 | end do |
---|
| 102 | end do |
---|
| 103 | end do |
---|
[72] | 104 | !$OMP END DO |
---|
[4] | 105 | |
---|
| 106 | ! *************************** Z VELOCITIES ****************** |
---|
| 107 | |
---|
[72] | 108 | !$OMP DO |
---|
[4] | 109 | do j=2,ny-1 |
---|
| 110 | do i=2,nx-1 |
---|
| 111 | |
---|
| 112 | divu(i,j)=((uxbar(i+1,j)*hmx(i+1,j)-uxbar(i,j)*hmx(i,j)) & ! formulation explicite |
---|
| 113 | + (uybar(i,j+1)*hmy(i,j+1)-uybar(i,j)*hmy(i,j)))/dx |
---|
| 114 | |
---|
| 115 | uzr(i,j,nz)=bmelt(i,j) |
---|
| 116 | xx(i,j)=uzr(i,j,1) |
---|
| 117 | hdd=bm(i,j)-bmelt(i,j)-divu(i,j) ! formulation explicite de dH/dt |
---|
| 118 | |
---|
| 119 | ! ------------------------------------------------------------------- |
---|
| 120 | ! attention uzr contient maintenant : |
---|
| 121 | ! uzr=uz +u ((1-e)*dH/dx+dB/dx)+ ((1-e)*dH/dt+dB/dt) |
---|
| 122 | |
---|
| 123 | do k=1,nz-1 |
---|
| 124 | if (.not.flot(i,j).and.(front(i,j).eq.4.).and..not.flgzmx(i,j) & |
---|
| 125 | .and..not.flgzmx(i+1,j).and..not.flgzmy(i,j) & |
---|
| 126 | .and..not.flgzmy(i,j+1)) then |
---|
| 127 | uzr(i,j,k)=uzr(i,j,nz)-((sux(i+1,j,k)-sux(i,j,k)) & |
---|
| 128 | +(suy(i,j+1,k)-suy(i,j,k)))/dx |
---|
| 129 | uzr(i,j,k)=uzr(i,j,k)+hdd*cde(k) |
---|
| 130 | |
---|
| 131 | else |
---|
| 132 | |
---|
| 133 | ! dans l'ice-shelf la vitesse uzr est lineaire ... |
---|
| 134 | ! ------------------------------------------------- |
---|
| 135 | |
---|
| 136 | ! attention, dans les zones stream, le vitesse verticale n'est pas |
---|
| 137 | ! lineaire si elle est voisine d'un point .not.flgzmx |
---|
| 138 | ! le calcul se fait alors par differences finies |
---|
| 139 | |
---|
| 140 | uzr(i,j,k)= bm(i,j)+(k-1.)/(nz-1.)*(bmelt(i,j)-bm(i,j)) |
---|
| 141 | end if |
---|
| 142 | |
---|
| 143 | end do |
---|
| 144 | |
---|
| 145 | end do |
---|
| 146 | end do |
---|
[72] | 147 | !$OMP END DO |
---|
[4] | 148 | |
---|
[72] | 149 | !$OMP DO |
---|
[4] | 150 | do j=2,ny-1 |
---|
| 151 | do i=2,nx-1 |
---|
| 152 | uzsdot(i,j)=(uzr(i,j,1)-xx(i,j))/dtt |
---|
| 153 | uzk(i,j)=(-bdot(i,j)-hdot(i,j)+bm(i,j)-bmelt(i,j)) |
---|
| 154 | end do |
---|
| 155 | end do |
---|
[72] | 156 | !$OMP END DO |
---|
| 157 | !$OMP END PARALLEL |
---|
[4] | 158 | |
---|
| 159 | end subroutine SIA_velocities |
---|