source: codes/icosagcm/trunk/src/etat0_williamson.f90 @ 195

Last change on this file since 195 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM
File size: 4.9 KB
Line 
1MODULE 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, etat0_williamson_new
10   
11CONTAINS
12 
13 
14  SUBROUTINE etat0_williamson(f_h,f_u) ! Deprecated / to be removed
15  USE icosa
16  IMPLICIT NONE
17    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 :: ind
23   
24    DO ind=1,ndomain
25      IF (.NOT. assigned_domain(ind)) CYCLE
26      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    ENDDO
32
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
53    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'
55       STOP
56    END IF
57
58    IF(llm>1) THEN
59       IF(is_mpi_root) PRINT *, 'etat0_type=williamson91.5 (Williamson,1991) must be used with llm=1 but llm =',llm
60       STOP
61    END IF
62
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(:,:)
74    ENDDO
75
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
90
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
171 
172END MODULE etat0_williamson_mod
Note: See TracBrowser for help on using the repository browser.