1 | MODULE etat0_williamson_mod |
---|
2 | USE icosa |
---|
3 | PRIVATE |
---|
4 | REAL(rstd), PARAMETER :: h0=8.E3 |
---|
5 | REAL(rstd), PARAMETER :: R0=4 |
---|
6 | REAL(rstd), PARAMETER :: K0=7.848E-6 |
---|
7 | REAL(rstd), PARAMETER :: omega0=K0 |
---|
8 | |
---|
9 | PUBLIC etat0_williamson, compute_etat0_williamson |
---|
10 | CONTAINS |
---|
11 | |
---|
12 | |
---|
13 | SUBROUTINE etat0_williamson(f_h,f_u) |
---|
14 | USE icosa |
---|
15 | IMPLICIT NONE |
---|
16 | TYPE(t_field),POINTER :: f_h(:) |
---|
17 | TYPE(t_field),POINTER :: f_u(:) |
---|
18 | |
---|
19 | REAL(rstd),POINTER :: h(:) |
---|
20 | REAL(rstd),POINTER :: u(:) |
---|
21 | INTEGER :: ind |
---|
22 | |
---|
23 | DO ind=1,ndomain |
---|
24 | CALL swap_dimensions(ind) |
---|
25 | CALL swap_geometry(ind) |
---|
26 | h=f_h(ind) |
---|
27 | u=f_u(ind) |
---|
28 | CALL compute_etat0_williamson(h, u) |
---|
29 | ENDDO |
---|
30 | |
---|
31 | END SUBROUTINE etat0_williamson |
---|
32 | |
---|
33 | SUBROUTINE compute_etat0_williamson(hi, ue) |
---|
34 | USE icosa |
---|
35 | IMPLICIT NONE |
---|
36 | REAL(rstd),INTENT(OUT) :: hi(iim*jjm) |
---|
37 | REAL(rstd),INTENT(OUT) :: ue(3*iim*jjm) |
---|
38 | REAL(rstd) :: lon, lat |
---|
39 | REAL(rstd) :: nx(3),n_norm,Velocity(3) |
---|
40 | REAL(rstd) :: A,B,C |
---|
41 | REAL(rstd) :: v1(3),v2(3),ny(3) |
---|
42 | REAL(rstd) :: de_min |
---|
43 | |
---|
44 | INTEGER :: i,j,n |
---|
45 | |
---|
46 | |
---|
47 | DO j=jj_begin-1,jj_end+1 |
---|
48 | DO i=ii_begin-1,ii_end+1 |
---|
49 | n=(j-1)*iim+i |
---|
50 | CALL xyz2lonlat(xyz_i(n,:),lon,lat) |
---|
51 | A= 0.5*omega0*(2*omega+omega0)*cos(lat)**2 & |
---|
52 | + 0.25*K0**2*cos(lat)**(2*R0)*((R0+1)*cos(lat)**2+(2*R0**2-R0-2)-2*R0**2/cos(lat)**2) |
---|
53 | B=2*(omega+omega0)*K0/((R0+1)*(R0+2))*cos(lat)**R0*((R0**2+2*R0+2)-(R0+1)**2*cos(lat)**2) |
---|
54 | C=0.25*K0**2*cos(lat)**(2*R0)*((R0+1)*cos(lat)**2-(R0+2)) |
---|
55 | |
---|
56 | hi(n)=(g*h0+radius**2*A+radius**2*B*cos(R0*lon)+radius**2*C*cos(2*R0*lon))/g |
---|
57 | |
---|
58 | |
---|
59 | CALL compute_velocity(xyz_e(n+u_right,:),velocity) |
---|
60 | CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),nx) |
---|
61 | |
---|
62 | ue(n+u_right)=1e-10 |
---|
63 | n_norm=sqrt(sum(nx(:)**2)) |
---|
64 | IF (n_norm>1e-30) THEN |
---|
65 | nx=-nx/n_norm*ne(n,right) |
---|
66 | ue(n+u_right)=sum(nx(:)*velocity(:)) |
---|
67 | IF (ABS(ue(n+u_right))<1e-100) PRINT *,"ue(n+u_right) ==0",i,j,velocity(:) |
---|
68 | ENDIF |
---|
69 | |
---|
70 | |
---|
71 | |
---|
72 | CALL compute_velocity(xyz_e(n+u_lup,:),velocity) |
---|
73 | CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),nx) |
---|
74 | |
---|
75 | ue(n+u_lup)=1e-10 |
---|
76 | n_norm=sqrt(sum(nx(:)**2)) |
---|
77 | IF (n_norm>1e-30) THEN |
---|
78 | nx=-nx/n_norm*ne(n,lup) |
---|
79 | ue(n+u_lup)=sum(nx(:)*velocity(:)) |
---|
80 | IF (ABS(ue(n+u_lup))<1e-100) PRINT *,"ue(n+u_lup) ==0",i,j,velocity(:) |
---|
81 | ENDIF |
---|
82 | |
---|
83 | |
---|
84 | CALL compute_velocity(xyz_e(n+u_ldown,:),velocity) |
---|
85 | CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),nx) |
---|
86 | |
---|
87 | ue(n+u_ldown)=1e-10 |
---|
88 | n_norm=sqrt(sum(nx(:)**2)) |
---|
89 | IF (n_norm>1e-30) THEN |
---|
90 | nx=-nx/n_norm*ne(n,ldown) |
---|
91 | ue(n+u_ldown)=sum(nx(:)*velocity(:)) |
---|
92 | IF (ABS(ue(n+u_ldown))<1e-100) PRINT *,"ue(n+u_ldown) ==0",i,j |
---|
93 | ENDIF |
---|
94 | |
---|
95 | |
---|
96 | ENDDO |
---|
97 | ENDDO |
---|
98 | |
---|
99 | de_min=1e10 |
---|
100 | DO j=jj_begin,jj_end |
---|
101 | DO i=ii_begin,ii_end |
---|
102 | n=(j-1)*iim+i |
---|
103 | de_min=MIN(de_min,de(n+u_right),de(n+u_rup),de(n+u_lup),de(n+u_left),de(n+u_ldown),de(n+u_rdown)) |
---|
104 | ENDDO |
---|
105 | ENDDO |
---|
106 | PRINT *,"-----> de min :",de_min |
---|
107 | |
---|
108 | CONTAINS |
---|
109 | |
---|
110 | SUBROUTINE compute_velocity(x,velocity) |
---|
111 | IMPLICIT NONE |
---|
112 | REAL(rstd),INTENT(IN) :: x(3) |
---|
113 | REAL(rstd),INTENT(OUT) :: velocity(3) |
---|
114 | REAL(rstd) :: e_lat(3), e_lon(3) |
---|
115 | REAL(rstd) :: lon,lat |
---|
116 | REAL(rstd) :: u,v |
---|
117 | |
---|
118 | CALL xyz2lonlat(x/radius,lon,lat) |
---|
119 | e_lat(1) = -cos(lon)*sin(lat) |
---|
120 | e_lat(2) = -sin(lon)*sin(lat) |
---|
121 | e_lat(3) = cos(lat) |
---|
122 | |
---|
123 | e_lon(1) = -sin(lon) |
---|
124 | e_lon(2) = cos(lon) |
---|
125 | e_lon(3) = 0 |
---|
126 | |
---|
127 | u=radius*omega0*cos(lat)+radius*K0*cos(lat)**(R0-1)*(R0*sin(lat)**2-cos(lat)**2)*cos(R0*lon) |
---|
128 | v=-radius*K0*R0*cos(lat)**(R0-1)*sin(lat)*sin(R0*lon) |
---|
129 | |
---|
130 | Velocity=(u*e_lon+v*e_lat+1e-50) |
---|
131 | |
---|
132 | END SUBROUTINE compute_velocity |
---|
133 | |
---|
134 | END SUBROUTINE compute_etat0_williamson |
---|
135 | |
---|
136 | END MODULE etat0_williamson_mod |
---|