[4] | 1 | !> \file taubed-0.3.f90 |
---|
| 2 | !! Calcul la charge en chaque point de la grille puis appel la routine litho pour calculer la contribution |
---|
| 3 | !! de chaque point a la deflexion de la lithosphere |
---|
| 4 | !< |
---|
| 5 | |
---|
| 6 | !> SUBROUTINE: taubed() |
---|
| 7 | !! \author ... |
---|
| 8 | !! \date 16 Novembre 1999 |
---|
| 9 | !! @note Routine qui calcul la charge en chaque point de la grille |
---|
| 10 | !! puis appel la routine litho pour calculer la contribution |
---|
| 11 | !! de chaque point a la deflexion de la lithosphere |
---|
| 12 | !! @note En entree |
---|
| 13 | !! - H, |
---|
| 14 | !! - BSOC, |
---|
| 15 | !! - SEALEVEL, |
---|
| 16 | !! - RO, |
---|
| 17 | !! - ROW |
---|
| 18 | !! @note En sortie |
---|
| 19 | !! - CHARGE(1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC) : poids par unite de surface |
---|
| 20 | !! (unite ?) Elle est calculee initialement dans initial2 |
---|
| 21 | !! Poids de la colonne d'eau ou de la colonne de glace. |
---|
| 22 | !! a l'exterieur du domaine : 1-LBLOC:1 et NX+1:NX+LBLOC |
---|
| 23 | !! on donne les valeurs des bords de la grille |
---|
| 24 | !! CHARGE est utilise par litho uniquement |
---|
| 25 | !! - W1(NX,NY) est l'enfoncement courant, c'est le resultat |
---|
| 26 | !! de la routine litho |
---|
| 27 | !! W1 peut etre calcule de plusieurs facons selon le modele |
---|
| 28 | !! d'isostasie utilise |
---|
| 29 | !! @note Used module: |
---|
| 30 | !! @note - use module3D_phy |
---|
| 31 | !! @note - use param_phy_mod |
---|
| 32 | !! @note - use ISO_DECLAR |
---|
| 33 | !! |
---|
| 34 | !< |
---|
| 35 | |
---|
| 36 | !! **************************************************************** |
---|
| 37 | !! * BEDROCK ADJUSTMENT avec temps de reaction * |
---|
| 38 | !! * changement de nom B -> Bsoc * |
---|
| 39 | !! **************************************************************** |
---|
| 40 | |
---|
| 41 | |
---|
| 42 | |
---|
| 43 | subroutine taubed() |
---|
| 44 | |
---|
[74] | 45 | !$USE OMP_LIB |
---|
[4] | 46 | USE module3D_phy |
---|
| 47 | USE param_phy_mod |
---|
| 48 | USE ISO_DECLAR ! module de declaration des variables de l'isostasie |
---|
| 49 | |
---|
| 50 | implicit none |
---|
| 51 | |
---|
| 52 | ! ********* calcul de W1 l'enfoncement d'equilibre au temps t |
---|
| 53 | ! NLITH est defini dans isostasie et permet le choix du modele d'isostasie |
---|
| 54 | |
---|
[74] | 55 | |
---|
[4] | 56 | if (NLITH.eq.1) then |
---|
| 57 | ! avec rigidite de la lithosphere |
---|
[74] | 58 | !$OMP PARALLEL |
---|
| 59 | !$OMP DO |
---|
[4] | 60 | do J=1,NY |
---|
| 61 | do I=1,NX |
---|
[237] | 62 | if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then |
---|
[4] | 63 | ! glace ou terre |
---|
| 64 | CHARGE(I,J)=ROG*H(I,J) |
---|
| 65 | else |
---|
| 66 | ! ocean |
---|
[237] | 67 | CHARGE(I,J)=-ROWG*(BSOC(I,J)-SEALEVEL_2D(i,j)) |
---|
[4] | 68 | endif |
---|
| 69 | end do |
---|
| 70 | end do |
---|
[74] | 71 | !$OMP END DO |
---|
[4] | 72 | |
---|
| 73 | ! il faut remplir CHARGE dans les parties a l'exterieur de la grille : |
---|
| 74 | ! a l'exterieur de la grille CHARGE est egale a la valeur sur le bord |
---|
[74] | 75 | !$OMP DO |
---|
[4] | 76 | do J=1,NY |
---|
| 77 | CHARGE(1-LBLOC:0,J)=CHARGE(1,J) ! bord W |
---|
| 78 | CHARGE(NX+1:NX+LBLOC,J)=CHARGE(NX,J) ! bord E |
---|
| 79 | end do |
---|
[74] | 80 | !$OMP END DO |
---|
| 81 | !$OMP DO |
---|
[4] | 82 | do I=1,NX |
---|
| 83 | CHARGE(I,1-LBLOC:0)=CHARGE(I,1) ! bord S |
---|
| 84 | CHARGE(I,NY+1:NY+LBLOC)=CHARGE(I,NY) ! bord N |
---|
| 85 | end do |
---|
[74] | 86 | !$OMP END DO |
---|
[4] | 87 | |
---|
| 88 | ! valeur dans les quatres coins exterieurs au domaine |
---|
[74] | 89 | !$OMP WORKSHARE |
---|
[4] | 90 | CHARGE(1-LBLOC:0,1-LBLOC:0)=CHARGE(1,1) ! coin SW |
---|
| 91 | CHARGE(1-LBLOC:0,NY+1:NY+LBLOC)=CHARGE(1,NY) ! coin NW |
---|
| 92 | CHARGE(NX+1:NX+LBLOC,1-LBLOC:0)=CHARGE(NX,1) ! coin SE |
---|
| 93 | CHARGE(NX+1:NX+LBLOC,NY+1:NY+LBLOC)=CHARGE(NX,NY) ! coin NE |
---|
[74] | 94 | !$OMP END WORKSHARE |
---|
| 95 | !$OMP END PARALLEL |
---|
[4] | 96 | call litho |
---|
| 97 | |
---|
| 98 | else |
---|
| 99 | ! enfoncement local |
---|
[74] | 100 | !$OMP PARALLEL |
---|
| 101 | !$OMP DO |
---|
[4] | 102 | do J=1,NY |
---|
| 103 | do I=1,NX |
---|
[237] | 104 | if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then |
---|
[4] | 105 | ! glace ou terre |
---|
| 106 | W1(I,J)=RO/ROM*H(I,J) |
---|
| 107 | else |
---|
| 108 | ! ocean |
---|
[237] | 109 | W1(I,J)=-ROW/ROM*(BSOC(I,J)-SEALEVEL_2D(i,j)) |
---|
[4] | 110 | endif |
---|
| 111 | end do |
---|
| 112 | end do |
---|
[74] | 113 | !$OMP END DO |
---|
| 114 | !$OMP END PARALLEL |
---|
[4] | 115 | endif |
---|
| 116 | |
---|
| 117 | ! decroissance exponentielle de l'enfoncement |
---|
[74] | 118 | !$OMP PARALLEL |
---|
| 119 | !$OMP DO |
---|
[4] | 120 | do J=1,NY |
---|
| 121 | do I=1,NX |
---|
| 122 | BDOT(I,J)=((Bsoc0(I,J)-BSOC(I,J))-(W1(I,J)-W0(I,J)))/TAUSOC |
---|
[181] | 123 | BSOC(I,J)=BSOC(I,J)+BDOT(I,J)*DT_ISO |
---|
[4] | 124 | end do |
---|
| 125 | end do |
---|
[74] | 126 | !$OMP END DO |
---|
| 127 | !$OMP END PARALLEL |
---|
[4] | 128 | |
---|
| 129 | end subroutine taubed |
---|