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