| 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 getin_etat0, compute_etat0 |
|---|
| 10 | |
|---|
| 11 | CONTAINS |
|---|
| 12 | |
|---|
| 13 | SUBROUTINE getin_etat0 |
|---|
| 14 | USE mpipara, ONLY : is_mpi_root |
|---|
| 15 | USE disvert_mod, ONLY : caldyn_eta, eta_lag |
|---|
| 16 | IF(caldyn_eta /= eta_lag) THEN |
|---|
| 17 | IF(is_mpi_root) PRINT *, 'etat0_type=williamson91.6 (Williamson,1991) must be used with caldyn_eta=eta_lag' |
|---|
| 18 | STOP |
|---|
| 19 | END IF |
|---|
| 20 | |
|---|
| 21 | IF(llm>1) THEN |
|---|
| 22 | IF(is_mpi_root) PRINT *, 'etat0_type=williamson91.6 (Williamson,1991) must be used with llm=1 but llm =',llm |
|---|
| 23 | STOP |
|---|
| 24 | END IF |
|---|
| 25 | END SUBROUTINE getin_etat0 |
|---|
| 26 | |
|---|
| 27 | SUBROUTINE compute_etat0(ngrid,lon,lat, phis,mass,rhodz,ulon,ulat) |
|---|
| 28 | IMPLICIT NONE |
|---|
| 29 | INTEGER, INTENT(IN) :: ngrid |
|---|
| 30 | REAL(rstd),INTENT(IN) :: lon(ngrid) |
|---|
| 31 | REAL(rstd),INTENT(IN) :: lat(ngrid) |
|---|
| 32 | REAL(rstd),INTENT(OUT) :: phis(ngrid) |
|---|
| 33 | REAL(rstd),INTENT(OUT) :: mass(ngrid) |
|---|
| 34 | REAL(rstd),INTENT(OUT) :: rhodz(ngrid) |
|---|
| 35 | REAL(rstd),INTENT(OUT) :: ulon(ngrid) |
|---|
| 36 | REAL(rstd),INTENT(OUT) :: ulat(ngrid) |
|---|
| 37 | |
|---|
| 38 | REAL(rstd) :: coslat,sinlat, A,B,C |
|---|
| 39 | INTEGER :: ij |
|---|
| 40 | |
|---|
| 41 | DO ij=1,ngrid |
|---|
| 42 | coslat=cos(lat(ij)) |
|---|
| 43 | sinlat=sin(lat(ij)) |
|---|
| 44 | A = 0.5*omega0*(2*omega+omega0)*coslat**2 & |
|---|
| 45 | + 0.25*K0**2*coslat**(2*R0)*((R0+1)*coslat**2+(2*R0**2-R0-2)-2*R0**2/coslat**2) |
|---|
| 46 | B = 2*(omega+omega0)*K0/((R0+1)*(R0+2))*coslat**R0*((R0**2+2*R0+2)-(R0+1)**2*coslat**2) |
|---|
| 47 | C = 0.25*K0**2*coslat**(2*R0)*((R0+1)*coslat**2-(R0+2)) |
|---|
| 48 | |
|---|
| 49 | mass(ij) = (g*h0+radius**2*A+radius**2*B*cos(R0*lon(ij))+radius**2*C*cos(2*R0*lon(ij)))/g |
|---|
| 50 | ulon(ij) = radius*omega0*coslat+radius*K0*coslat**(R0-1)*(R0*sinlat**2-coslat**2)*cos(R0*lon(ij)) |
|---|
| 51 | ulat(ij) = -radius*K0*R0*coslat**(R0-1)*sinlat*sin(R0*lon(ij)) |
|---|
| 52 | ENDDO |
|---|
| 53 | |
|---|
| 54 | phis(:)=0. |
|---|
| 55 | rhodz(:)=mass(:) |
|---|
| 56 | |
|---|
| 57 | END SUBROUTINE compute_etat0 |
|---|
| 58 | |
|---|
| 59 | END MODULE etat0_williamson_mod |
|---|
