![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 204 for codes/icosagcm/trunk/src/etat0_williamson.f90
- Timestamp:
- 07/09/14 09:36:15 (10 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/trunk/src/etat0_williamson.f90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/etat0_williamson.f90
r186 r204 2 2 USE icosa 3 3 PRIVATE 4 5 6 7 8 9 PUBLIC etat0_williamson, etat0_williamson_new10 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 11 CONTAINS 12 13 14 SUBROUTINE etat0_williamson(f_h,f_u) ! Deprecated / to be removed15 USE icosa16 IMPLICIT NONE17 TYPE(t_field),POINTER :: f_h(:)18 TYPE(t_field),POINTER :: f_u(:)19 20 REAL(rstd),POINTER :: h(:)21 REAL(rstd),POINTER :: u(:)22 INTEGER :: ind23 24 DO ind=1,ndomain25 IF (.NOT. assigned_domain(ind)) CYCLE26 CALL swap_dimensions(ind)27 CALL swap_geometry(ind)28 h=f_h(ind)29 u=f_u(ind)30 CALL compute_etat0_williamson(h, u)31 ENDDO32 12 33 END SUBROUTINE etat0_williamson 34 35 SUBROUTINE etat0_williamson_new(f_phis,f_mass,f_theta_rhodz,f_u, f_q) 36 USE icosa 37 USE mpipara, ONLY : is_mpi_root 38 USE disvert_mod, ONLY : caldyn_eta, eta_lag 39 40 IMPLICIT NONE 41 TYPE(t_field),POINTER :: f_phis(:) 42 TYPE(t_field),POINTER :: f_mass(:) 43 TYPE(t_field),POINTER :: f_theta_rhodz(:) 44 TYPE(t_field),POINTER :: f_u(:) 45 TYPE(t_field),POINTER :: f_q(:) 46 47 REAL(rstd),POINTER :: phis(:) 48 REAL(rstd),POINTER :: h(:,:) 49 REAL(rstd),POINTER :: theta_rhodz(:,:) 50 REAL(rstd),POINTER :: u(:,:) 51 INTEGER :: ind 52 13 SUBROUTINE getin_etat0 14 USE mpipara, ONLY : is_mpi_root 15 USE disvert_mod, ONLY : caldyn_eta, eta_lag 53 16 IF(caldyn_eta /= eta_lag) THEN 54 IF(is_mpi_root) PRINT *, 'etat0_type=williamson91. 5(Williamson,1991) must be used with caldyn_eta=eta_lag'17 IF(is_mpi_root) PRINT *, 'etat0_type=williamson91.6 (Williamson,1991) must be used with caldyn_eta=eta_lag' 55 18 STOP 56 19 END IF 57 20 58 21 IF(llm>1) THEN 59 IF(is_mpi_root) PRINT *, 'etat0_type=williamson91. 5(Williamson,1991) must be used with llm=1 but llm =',llm22 IF(is_mpi_root) PRINT *, 'etat0_type=williamson91.6 (Williamson,1991) must be used with llm=1 but llm =',llm 60 23 STOP 61 24 END IF 25 END SUBROUTINE getin_etat0 62 26 63 DO ind=1,ndomain 64 IF (.NOT. assigned_domain(ind)) CYCLE 65 CALL swap_dimensions(ind) 66 CALL swap_geometry(ind) 67 h=f_mass(ind) 68 u=f_u(ind) 69 theta_rhodz=f_theta_rhodz(ind) 70 phis=f_phis(ind) 71 CALL compute_etat0_williamson(h(:,1), u(:,1)) 72 phis(:)=0. 73 theta_rhodz(:,:) = h(:,:) 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)) 74 52 ENDDO 75 53 76 END SUBROUTINE etat0_williamson_new 77 78 SUBROUTINE compute_etat0_williamson(hi, ue) 79 USE icosa 80 IMPLICIT NONE 81 REAL(rstd),INTENT(OUT) :: hi(iim*jjm) 82 REAL(rstd),INTENT(OUT) :: ue(3*iim*jjm) 83 REAL(rstd) :: lon, lat 84 REAL(rstd) :: nx(3),n_norm,Velocity(3) 85 REAL(rstd) :: A,B,C 86 REAL(rstd) :: v1(3),v2(3),ny(3) 87 REAL(rstd) :: de_min 88 89 INTEGER :: i,j,n 54 phis(:)=0. 55 rhodz(:)=mass(:) 90 56 91 92 DO j=jj_begin-1,jj_end+1 93 DO i=ii_begin-1,ii_end+1 94 n=(j-1)*iim+i 95 CALL xyz2lonlat(xyz_i(n,:),lon,lat) 96 A= 0.5*omega0*(2*omega+omega0)*cos(lat)**2 & 97 + 0.25*K0**2*cos(lat)**(2*R0)*((R0+1)*cos(lat)**2+(2*R0**2-R0-2)-2*R0**2/cos(lat)**2) 98 B=2*(omega+omega0)*K0/((R0+1)*(R0+2))*cos(lat)**R0*((R0**2+2*R0+2)-(R0+1)**2*cos(lat)**2) 99 C=0.25*K0**2*cos(lat)**(2*R0)*((R0+1)*cos(lat)**2-(R0+2)) 100 101 hi(n)=(g*h0+radius**2*A+radius**2*B*cos(R0*lon)+radius**2*C*cos(2*R0*lon))/g 102 103 104 CALL compute_velocity(xyz_e(n+u_right,:),velocity) 105 CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),nx) 106 107 ue(n+u_right)=1e-10 108 n_norm=sqrt(sum(nx(:)**2)) 109 IF (n_norm>1e-30) THEN 110 nx=-nx/n_norm*ne(n,right) 111 ue(n+u_right)=sum(nx(:)*velocity(:)) 112 IF (ABS(ue(n+u_right))<1e-100) PRINT *,"ue(n+u_right) ==0",i,j,velocity(:) 113 ENDIF 114 115 116 117 CALL compute_velocity(xyz_e(n+u_lup,:),velocity) 118 CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),nx) 119 120 ue(n+u_lup)=1e-10 121 n_norm=sqrt(sum(nx(:)**2)) 122 IF (n_norm>1e-30) THEN 123 nx=-nx/n_norm*ne(n,lup) 124 ue(n+u_lup)=sum(nx(:)*velocity(:)) 125 IF (ABS(ue(n+u_lup))<1e-100) PRINT *,"ue(n+u_lup) ==0",i,j,velocity(:) 126 ENDIF 127 128 129 CALL compute_velocity(xyz_e(n+u_ldown,:),velocity) 130 CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),nx) 131 132 ue(n+u_ldown)=1e-10 133 n_norm=sqrt(sum(nx(:)**2)) 134 IF (n_norm>1e-30) THEN 135 nx=-nx/n_norm*ne(n,ldown) 136 ue(n+u_ldown)=sum(nx(:)*velocity(:)) 137 IF (ABS(ue(n+u_ldown))<1e-100) PRINT *,"ue(n+u_ldown) ==0",i,j 138 ENDIF 139 140 141 ENDDO 142 ENDDO 143 144 CONTAINS 145 146 SUBROUTINE compute_velocity(x,velocity) 147 IMPLICIT NONE 148 REAL(rstd),INTENT(IN) :: x(3) 149 REAL(rstd),INTENT(OUT) :: velocity(3) 150 REAL(rstd) :: e_lat(3), e_lon(3) 151 REAL(rstd) :: lon,lat 152 REAL(rstd) :: u,v 153 154 CALL xyz2lonlat(x/radius,lon,lat) 155 e_lat(1) = -cos(lon)*sin(lat) 156 e_lat(2) = -sin(lon)*sin(lat) 157 e_lat(3) = cos(lat) 158 159 e_lon(1) = -sin(lon) 160 e_lon(2) = cos(lon) 161 e_lon(3) = 0 162 163 u=radius*omega0*cos(lat)+radius*K0*cos(lat)**(R0-1)*(R0*sin(lat)**2-cos(lat)**2)*cos(R0*lon) 164 v=-radius*K0*R0*cos(lat)**(R0-1)*sin(lat)*sin(R0*lon) 165 166 Velocity=(u*e_lon+v*e_lat+1e-50) 167 168 END SUBROUTINE compute_velocity 169 170 END SUBROUTINE compute_etat0_williamson 57 END SUBROUTINE compute_etat0 171 58 172 59 END MODULE etat0_williamson_mod
Note: See TracChangeset
for help on using the changeset viewer.
