[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 | ! ========================================================================== |
---|
| 25 | |
---|
| 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 | |
---|
| 34 | |
---|
| 35 | do k=1,nz |
---|
| 36 | do j=2,ny-1 |
---|
| 37 | do i=2,nx-1 |
---|
| 38 | |
---|
| 39 | ! flgzmx pour tous les points qui sont traites par l'equation elliptique |
---|
| 40 | ! diagnoshelf |
---|
| 41 | |
---|
| 42 | ! Selon X |
---|
| 43 | ! ------------------------------------- |
---|
| 44 | ! if (.not.flgzmx(i,j)) then ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) |
---|
| 45 | ! ux(i,j,k)=ddbx(i,j) |
---|
| 46 | ! sux(i,j,k)=ddbx(i,j)*cde(k) |
---|
| 47 | |
---|
| 48 | if (.not.flotmx(i,j)) then ! version aout 2010 (Cat) |
---|
| 49 | ! utilise ubx au lieu de ddbx (compatibilite avec mix et spinup) |
---|
| 50 | ! et agit sur tous les points poses |
---|
| 51 | |
---|
| 52 | ux(i,j,k) = 0. |
---|
| 53 | sux(i,j,k)= 0. |
---|
| 54 | |
---|
| 55 | do iglen=n1poly,n2poly |
---|
| 56 | |
---|
| 57 | ux(i,j,k) = ux(i,j,k)+ddx(i,j,iglen)*sa_mx(i,j,k,iglen) |
---|
| 58 | sux(i,j,k) = sux(i,j,k)+ddx(i,j,iglen)*s2a_mx(i,j,k,iglen) |
---|
| 59 | end do |
---|
| 60 | ux(i,j,k) = ux(i,j,k)*(-sdx(i,j)) |
---|
| 61 | ux(i,j,k) = ux(i,j,k) + ubx(i,j) ! ajoute le glissement |
---|
| 62 | |
---|
| 63 | sux(i,j,k) = sux(i,j,k)*sdx(i,j)*hmx(i,j) |
---|
| 64 | sux(i,j,k) = sux(i,j,k) - ubx(i,j) * cde(k) |
---|
| 65 | |
---|
| 66 | else |
---|
| 67 | sux(i,j,k) = (nz-k)*1./(nz-1.)*uxbar(i,j)*hmx(i,j) |
---|
| 68 | ux(i,j,k) = uxbar(i,j) |
---|
| 69 | end if |
---|
| 70 | |
---|
| 71 | ! Selon Y |
---|
| 72 | ! ------------------------------------- |
---|
| 73 | ! if (.not.flgzmy(i,j)) then |
---|
| 74 | ! uy(i,j,k)=ddby(i,j) |
---|
| 75 | ! suy(i,j,k)=ddby(i,j)*cde(k) |
---|
| 76 | |
---|
| 77 | if (.not.flotmy(i,j)) then ! version aout 2010 (Cat) |
---|
| 78 | ! utilise uby au lieu de ddby (compatibilite avec mix et spinup) |
---|
| 79 | ! et agit sur tous les points poses |
---|
| 80 | |
---|
| 81 | uy(i,j,k) = 0. |
---|
| 82 | suy(i,j,k)= 0. |
---|
| 83 | |
---|
| 84 | |
---|
| 85 | do iglen=n1poly,n2poly |
---|
| 86 | uy(i,j,k)=uy(i,j,k)+ddy(i,j,iglen)*sa_my(i,j,k,iglen) |
---|
| 87 | suy(i,j,k)=suy(i,j,k)+ddy(i,j,iglen)*s2a_my(i,j,k,iglen) |
---|
| 88 | end do |
---|
| 89 | uy(i,j,k) = uy(i,j,k)*(-sdy(i,j)) |
---|
| 90 | uy(i,j,k) = uy(i,j,k) + uby(i,j) ! ajoute le glissement |
---|
| 91 | |
---|
| 92 | suy(i,j,k) = suy(i,j,k)*sdy(i,j)*hmy(i,j) |
---|
| 93 | suy(i,j,k) = suy(i,j,k) - uby(i,j) * cde(k) |
---|
| 94 | |
---|
| 95 | else |
---|
| 96 | suy(i,j,k) = (nz-k)*1./(nz-1)*uybar(i,j)*hmy(i,j) |
---|
| 97 | uy(i,j,k)=uybar(i,j) |
---|
| 98 | endif |
---|
| 99 | |
---|
| 100 | end do |
---|
| 101 | end do |
---|
| 102 | end do |
---|
| 103 | |
---|
| 104 | |
---|
| 105 | ! *************************** Z VELOCITIES ****************** |
---|
| 106 | |
---|
| 107 | |
---|
| 108 | do j=2,ny-1 |
---|
| 109 | do i=2,nx-1 |
---|
| 110 | |
---|
| 111 | divu(i,j)=((uxbar(i+1,j)*hmx(i+1,j)-uxbar(i,j)*hmx(i,j)) & ! formulation explicite |
---|
| 112 | + (uybar(i,j+1)*hmy(i,j+1)-uybar(i,j)*hmy(i,j)))/dx |
---|
| 113 | |
---|
| 114 | uzr(i,j,nz)=bmelt(i,j) |
---|
| 115 | xx(i,j)=uzr(i,j,1) |
---|
| 116 | hdd=bm(i,j)-bmelt(i,j)-divu(i,j) ! formulation explicite de dH/dt |
---|
| 117 | |
---|
| 118 | ! ------------------------------------------------------------------- |
---|
| 119 | ! attention uzr contient maintenant : |
---|
| 120 | ! uzr=uz +u ((1-e)*dH/dx+dB/dx)+ ((1-e)*dH/dt+dB/dt) |
---|
| 121 | |
---|
| 122 | do k=1,nz-1 |
---|
| 123 | if (.not.flot(i,j).and.(front(i,j).eq.4.).and..not.flgzmx(i,j) & |
---|
| 124 | .and..not.flgzmx(i+1,j).and..not.flgzmy(i,j) & |
---|
| 125 | .and..not.flgzmy(i,j+1)) then |
---|
| 126 | uzr(i,j,k)=uzr(i,j,nz)-((sux(i+1,j,k)-sux(i,j,k)) & |
---|
| 127 | +(suy(i,j+1,k)-suy(i,j,k)))/dx |
---|
| 128 | uzr(i,j,k)=uzr(i,j,k)+hdd*cde(k) |
---|
| 129 | |
---|
| 130 | else |
---|
| 131 | |
---|
| 132 | ! dans l'ice-shelf la vitesse uzr est lineaire ... |
---|
| 133 | ! ------------------------------------------------- |
---|
| 134 | |
---|
| 135 | ! attention, dans les zones stream, le vitesse verticale n'est pas |
---|
| 136 | ! lineaire si elle est voisine d'un point .not.flgzmx |
---|
| 137 | ! le calcul se fait alors par differences finies |
---|
| 138 | |
---|
| 139 | uzr(i,j,k)= bm(i,j)+(k-1.)/(nz-1.)*(bmelt(i,j)-bm(i,j)) |
---|
| 140 | end if |
---|
| 141 | |
---|
| 142 | end do |
---|
| 143 | |
---|
| 144 | end do |
---|
| 145 | end do |
---|
| 146 | |
---|
| 147 | do j=2,ny-1 |
---|
| 148 | do i=2,nx-1 |
---|
| 149 | uzsdot(i,j)=(uzr(i,j,1)-xx(i,j))/dtt |
---|
| 150 | uzk(i,j)=(-bdot(i,j)-hdot(i,j)+bm(i,j)-bmelt(i,j)) |
---|
| 151 | end do |
---|
| 152 | end do |
---|
| 153 | |
---|
| 154 | |
---|
| 155 | end subroutine SIA_velocities |
---|