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 | |
---|
45 | !$USE OMP_LIB |
---|
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 | |
---|
55 | |
---|
56 | if (NLITH.eq.1) then |
---|
57 | ! avec rigidite de la lithosphere |
---|
58 | !$OMP PARALLEL |
---|
59 | !$OMP DO |
---|
60 | do J=1,NY |
---|
61 | do I=1,NX |
---|
62 | if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then |
---|
63 | ! glace ou terre |
---|
64 | CHARGE(I,J)=ROG*H(I,J) |
---|
65 | else |
---|
66 | ! ocean |
---|
67 | CHARGE(I,J)=-ROWG*(BSOC(I,J)-SEALEVEL_2D(i,j)) |
---|
68 | endif |
---|
69 | end do |
---|
70 | end do |
---|
71 | !$OMP END DO |
---|
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 |
---|
75 | !$OMP DO |
---|
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 |
---|
80 | !$OMP END DO |
---|
81 | !$OMP DO |
---|
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 |
---|
86 | !$OMP END DO |
---|
87 | |
---|
88 | ! valeur dans les quatres coins exterieurs au domaine |
---|
89 | !$OMP WORKSHARE |
---|
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 |
---|
94 | !$OMP END WORKSHARE |
---|
95 | !$OMP END PARALLEL |
---|
96 | call litho |
---|
97 | |
---|
98 | else |
---|
99 | ! enfoncement local |
---|
100 | !$OMP PARALLEL |
---|
101 | !$OMP DO |
---|
102 | do J=1,NY |
---|
103 | do I=1,NX |
---|
104 | if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then |
---|
105 | ! glace ou terre |
---|
106 | W1(I,J)=RO/ROM*H(I,J) |
---|
107 | else |
---|
108 | ! ocean |
---|
109 | W1(I,J)=-ROW/ROM*(BSOC(I,J)-SEALEVEL_2D(i,j)) |
---|
110 | endif |
---|
111 | end do |
---|
112 | end do |
---|
113 | !$OMP END DO |
---|
114 | !$OMP END PARALLEL |
---|
115 | endif |
---|
116 | |
---|
117 | ! decroissance exponentielle de l'enfoncement |
---|
118 | !$OMP PARALLEL |
---|
119 | !$OMP DO |
---|
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 |
---|
123 | BSOC(I,J)=BSOC(I,J)+BDOT(I,J)*DT_ISO |
---|
124 | end do |
---|
125 | end do |
---|
126 | !$OMP END DO |
---|
127 | !$OMP END PARALLEL |
---|
128 | |
---|
129 | end subroutine taubed |
---|