source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/physics.f90 @ 316

Last change on this file since 316 was 221, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM
File size: 6.5 KB
Line 
1MODULE physics_mod
2
3  USE field_mod
4
5  PRIVATE
6
7  INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2, phys_lmdz_generic=3
8
9  INTEGER :: nb_extra_physics_2D, nb_extra_physics_3D, phys_type
10  TYPE(t_field),POINTER :: f_extra_physics_2D(:), f_extra_physics_3D(:)
11
12  CHARACTER(LEN=255) :: physics_type="automatic"
13!$OMP THREADPRIVATE(physics_type)
14
15  PUBLIC :: physics, init_physics
16CONTAINS
17
18  SUBROUTINE init_physics
19    USE mpipara
20    USE etat0_mod
21    USE icosa
22    USE physics_interface_mod
23    USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
24    USE physics_lmdz_generic_mod, ONLY : init_physics_lmdz_generic=>init_physics
25    IMPLICIT NONE
26
27    physics_type='automatic'
28    CALL getin("physics",physics_type)
29
30    SELECT CASE(TRIM(physics_type))
31    CASE ('automatic')
32       etat0_type='jablonowsky06'
33       CALL getin("etat0",etat0_type)
34       SELECT CASE(TRIM(etat0_type))
35       CASE('held_suarez')
36          phys_type = phys_HS94
37       CASE DEFAULT
38          IF(is_mpi_root) PRINT*,"NO PHYSICAL PACKAGE USED"
39          phys_type = phys_none
40       END SELECT
41    CASE ('phys_lmdz_generic')
42       CALL init_physics_lmdz_generic
43       phys_type=phys_lmdz_generic
44    CASE ('dcmip')
45       CALL init_physics_dcmip
46       nb_extra_physics_2D=1 ! precl
47       nb_extra_physics_3D=0
48       phys_type = phys_DCMIP
49    CASE DEFAULT
50       IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
51            TRIM(physics_type), '> options are <automatic>, <dcmip>, <dry>'
52       STOP
53    END SELECT
54
55    IF(is_mpi_root) THEN
56       PRINT *, 'phys_type = ',phys_type
57       PRINT *, 'nb_extra_physics_2D = ', nb_extra_physics_2D
58       PRINT *, 'nb_extra_physics_3D = ', nb_extra_physics_3D
59    END IF
60    IF(nb_extra_physics_2D>0) CALL allocate_field(f_extra_physics_2D,field_t,type_real,nb_extra_physics_2D)
61    IF(nb_extra_physics_3D>0) CALL allocate_field(f_extra_physics_3D,field_t,type_real,llm,nb_extra_physics_3D)
62
63  END SUBROUTINE init_physics
64
65  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
66    USE icosa
67    USE physics_interface_mod
68    USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics
69 !   USE etat0_mod
70    USE etat0_heldsz_mod
71    IMPLICIT NONE
72    INTEGER, INTENT(IN)   :: it
73    TYPE(t_field),POINTER :: f_phis(:)
74    TYPE(t_field),POINTER :: f_ps(:)
75    TYPE(t_field),POINTER :: f_theta_rhodz(:)
76    TYPE(t_field),POINTER :: f_ue(:)
77    TYPE(t_field),POINTER :: f_wflux(:)
78    TYPE(t_field),POINTER :: f_q(:)
79    REAL(rstd),POINTER :: phis(:)
80    REAL(rstd),POINTER :: ps(:)
81    REAL(rstd),POINTER :: theta_rhodz(:,:)
82    REAL(rstd),POINTER :: ue(:,:)
83    REAL(rstd),POINTER :: q(:,:,:)
84
85    LOGICAL:: firstcall,lastcall
86    INTEGER :: ind
87    TYPE(physics_inout) :: args
88
89    SELECT CASE(phys_type)
90    CASE (phys_none)
91       ! No physics, do nothing
92    CASE(phys_HS94)
93       CALL held_suarez(f_ps,f_theta_rhodz,f_ue)
94    CASE (phys_lmdz_generic)
95       CALL physics_lmdz_generic(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
96    CASE DEFAULT
97       CALL transfert_request(f_ps,req_i1)
98       CALL transfert_request(f_theta_rhodz,req_i1)
99       CALL transfert_request(f_ue,req_e1_vect)
100       CALL transfert_request(f_q,req_i1)
101       
102       args%dt_phys = dt
103       DO ind=1,ndomain
104          IF (.NOT. assigned_domain(ind)) CYCLE
105          CALL swap_dimensions(ind)
106          CALL swap_geometry(ind)
107     
108          phis=f_phis(ind)
109          ps=f_ps(ind)
110          theta_rhodz=f_theta_rhodz(ind)
111          ue=f_ue(ind)
112          q=f_q(ind)
113
114          IF(nb_extra_physics_2D>0) args%extra_2D=f_extra_physics_2D(ind)
115          IF(nb_extra_physics_3D>0) args%extra_3D=f_extra_physics_3D(ind)
116          CALL physics_column(args, phis, ps, theta_rhodz, ue, q)
117       ENDDO
118
119       IF (mod(it,itau_out)==0 ) THEN
120          IF(nb_extra_physics_2D>0) CALL writefield("extra_physics_2D",f_extra_physics_2D)
121          IF(nb_extra_physics_3D>0) CALL writefield("extra_physics_3D",f_extra_physics_3D)
122       ENDIF
123    END SELECT
124
125  END SUBROUTINE physics
126
127  SUBROUTINE physics_column(args, phis, ps, theta_rhodz, ue, q)
128    USE icosa
129    USE wind_mod
130    USE pression_mod
131    USE theta2theta_rhodz_mod
132    USE physics_interface_mod
133    USE physics_dcmip_mod
134    IMPLICIT NONE
135    TYPE(physics_inout) :: args   
136    REAL(rstd) :: phis(iim*jjm)
137    REAL(rstd) :: ps(iim*jjm)
138    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
139    REAL(rstd) :: ue(3*iim*jjm,llm)
140    REAL(rstd), TARGET :: q(iim*jjm,llm,nqtot)
141    ! local arrays
142    REAL(rstd), TARGET :: lat(iim*jjm)
143    REAL(rstd), TARGET :: lon(iim*jjm)
144    REAL(rstd), TARGET :: p(iim*jjm,llm+1)
145    REAL(rstd), TARGET :: Temp(iim*jjm,llm)
146    REAL(rstd), TARGET :: ulon(iim*jjm,llm)
147    REAL(rstd), TARGET :: ulat(iim*jjm,llm)
148    REAL(rstd), TARGET :: dTemp(iim*jjm,llm)
149    REAL(rstd), TARGET :: dulon(iim*jjm,llm)
150    REAL(rstd), TARGET :: dulat(iim*jjm,llm)
151    REAL(rstd), TARGET :: dq(iim*jjm,llm,nqtot)
152    REAL(rstd) :: uc(iim*jjm,3,llm)  ! 3D velocity at cell centers
153
154    INTEGER :: i,j,ij,l
155    REAL(rstd) :: due, dt2
156
157    DO j=jj_begin,jj_end
158      DO i=ii_begin,ii_end
159        ij=(j-1)*iim+i
160        CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij)) 
161      ENDDO
162    ENDDO
163
164    ! Reconstruct wind vector at hexagons
165    CALL compute_pression(ps,p,0)
166    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,Temp,0)
167    CALL compute_wind_centered(ue,uc)
168    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
169    args%ngrid = iim*jjm
170    args%lon => lon
171    args%lat => lat
172    args%p => p
173    args%Temp => Temp
174    args%ulon => ulon
175    args%ulat => ulat
176    args%q => q
177    args%dTemp => dTemp
178    args%dulon => dulon
179    args%dulat => dulat
180    args%dq => dq
181
182    SELECT CASE(phys_type)
183      CASE (phys_DCMIP)
184         CALL compute_phys_wrap(args)
185    END SELECT
186
187    q = q + args%dt_phys * dq
188    Temp = Temp + args%dt_phys * dTemp
189    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
190   
191    ! Reconstruct wind tendencies at edges and add
192    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,uc)
193    dt2=.5*args%dt_phys
194    DO l=1,llm
195      DO j=jj_begin,jj_end
196        DO i=ii_begin,ii_end
197          ij=(j-1)*iim+i
198          due = sum((uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
199          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
200
201          due = sum( (uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
202          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
203
204          due = sum( (uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
205          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
206        ENDDO
207      ENDDO
208    ENDDO
209
210  END SUBROUTINE physics_column
211
212END MODULE physics_mod
Note: See TracBrowser for help on using the repository browser.