source: codes/icosagcm/trunk/src/physics_dcmip2016.f90 @ 381

Last change on this file since 381 was 381, checked in by ymipsl, 8 years ago

Add DCMIP2016 physics

(first guess)

YM

File size: 6.1 KB
Line 
1MODULE physics_dcmip2016_mod
2  USE ICOSA
3  PRIVATE
4 
5  INTEGER,SAVE :: testcase
6!$OMP THREADPRIVATE(testcase)
7
8  TYPE(t_field),POINTER :: f_out_i(:)
9  REAL(rstd),POINTER :: out_i(:,:)
10
11  TYPE(t_field),POINTER  :: f_precl(:)
12  REAL(rstd),ALLOCATABLE :: precl_packed(:)
13
14  PUBLIC :: init_physics, full_physics, write_physics
15
16  INTEGER, PARAMETER :: dry_baroclinic=0
17  INTEGER, PARAMETER :: moist_baroclinic=1
18  INTEGER, PARAMETER :: cyclone=2
19  INTEGER, PARAMETER :: supercell=3
20
21 
22CONTAINS
23
24  SUBROUTINE init_physics
25    USE physics_interface_mod
26    IMPLICIT NONE
27    INTEGER :: ngrid
28    CHARACTER(LEN=255) :: testcase_str 
29   
30    CALL getin("physics_dcmip2016",testcase_str)
31   
32    SELECT CASE (TRIM(testcase_str))
33      CASE ('dry_baroclinic')
34        testcase=dry_baroclinic
35      CASE ('moist_baroclinic') 
36        testcase=moist_baroclinic
37      CASE ('cyclone') 
38        testcase=cyclone
39      CASE ('supercell') 
40        testcase=supercell
41      CASE DEFAULT
42         PRINT*, 'Bad selector for dcmip2016 test case <', testcase_str, &
43              '> options are <dry_baroclinic>, <moist_baroclinic>, <cyclone>, <supercell>'
44         STOP
45    END SELECT
46
47
48    ngrid = physics_inout%ngrid
49    ! Input
50    ALLOCATE(physics_inout%Ai(ngrid))
51    ALLOCATE(physics_inout%lon(ngrid))
52    ALLOCATE(physics_inout%lat(ngrid))
53    ALLOCATE(physics_inout%phis(ngrid))
54    ALLOCATE(physics_inout%p(ngrid,llm+1))
55    ALLOCATE(physics_inout%pk(ngrid,llm))
56    ALLOCATE(physics_inout%Temp(ngrid,llm))
57    ALLOCATE(physics_inout%ulon(ngrid,llm))
58    ALLOCATE(physics_inout%ulat(ngrid,llm))
59    ALLOCATE(physics_inout%q(ngrid,llm,nqtot))
60    ! Output (tendencies)
61    ALLOCATE(physics_inout%dTemp(ngrid,llm))
62    ALLOCATE(physics_inout%dulon(ngrid,llm))
63    ALLOCATE(physics_inout%dulat(ngrid,llm))
64    ALLOCATE(physics_inout%dq(ngrid,llm,nqtot))
65    ! Physics-specific data
66    ALLOCATE(precl_packed(ngrid))
67    CALL allocate_field(f_precl, field_t,type_real)
68
69    PRINT *, 'init_physics_new', SIZE(physics_inout%Ai)
70  END SUBROUTINE init_physics
71
72  SUBROUTINE full_physics
73    USE physics_interface_mod
74    CALL compute_physics(physics_inout%ngrid, physics_inout%dt_phys, &
75         physics_inout%lon, physics_inout%lat, physics_inout%p, physics_inout%pk, physics_inout%Temp, & 
76         physics_inout%ulon, physics_inout%ulat, physics_inout%q(:,:,1:5), &
77         physics_inout%dTemp, physics_inout%dulon, physics_inout%dulat, &
78         physics_inout%dq(:,:,1:5), precl_packed)
79  END SUBROUTINE full_physics
80
81  SUBROUTINE write_physics
82    USE output_field_mod
83    USE physics_interface_mod
84    CALL unpack_field(f_precl, precl_packed)
85    CALL output_field("precl",f_precl)
86  END SUBROUTINE write_physics
87
88  SUBROUTINE compute_physics(ngrid,dt_phys,lon, lat, p, pk, Temp,u,v,q, dTemp,du,dv,dq, precl)
89    USE icosa
90    USE dcmip2016_simple_physics_mod
91    USE dcmip2016_kessler_physic_mod
92    USE terminator
93    IMPLICIT NONE
94    INTEGER    :: ngrid
95    REAL(rstd) :: lat(ngrid)
96    REAL(rstd) :: lon(ngrid)
97    REAL(rstd) :: ps(ngrid)
98    REAL(rstd) :: precl(ngrid)
99    ! arguments with bottom-up indexing (DYNAMICO)
100    REAL(rstd) :: p(ngrid,llm+1)
101    REAL(rstd) :: pk(ngrid,llm)
102    REAL(rstd) :: Temp(ngrid,llm)
103    REAL(rstd) :: u(ngrid,llm)
104    REAL(rstd) :: v(ngrid,llm)
105    REAL(rstd) :: q(ngrid,llm,5)
106    REAL(rstd) :: dTemp(ngrid,llm)
107    REAL(rstd) :: du(ngrid,llm)
108    REAL(rstd) :: dv(ngrid,llm)
109    REAL(rstd) :: dq(ngrid,llm,5)
110    ! local arrays with top-down vertical indexing (DCMIP)
111    REAL(rstd) :: pint(ngrid,llm+1)
112    REAL(rstd) :: pmid(ngrid,llm)
113    REAL(rstd) :: pdel(ngrid,llm)
114    REAL(rstd) :: Tfi(ngrid,llm)
115    REAL(rstd) :: ufi(ngrid,llm)
116    REAL(rstd) :: vfi(ngrid,llm)
117    REAL(rstd) :: qfi(ngrid,llm,5)
118
119    REAL(rstd)  :: rho(llm), z(llm), theta(llm)
120    REAL(rstd)  :: lastz
121    REAL(rstd)  :: dcl1,dcl2
122    INTEGER :: l,ll,ij
123    REAL(rstd) :: dt_phys, inv_dt
124
125    ! prepare input fields and mirror vertical index     
126    ps(:) = p(:,1) ! surface pressure
127
128    DO l=1,llm+1
129      DO ij=1,ngrid
130          pint(ij,l)=p(ij,llm+2-l)
131      ENDDO
132    ENDDO
133
134    DO l=1,llm
135       ll=llm+1-l
136       DO ij=1,ngrid
137          pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) ! Pressure inside layers
138          pdel(ij,l)=pint(ij,l+1)-pint(ij,l)       ! Pressure difference between two layers
139          ufi(ij,l)=u(ij,ll)
140          vfi(ij,l)=v(ij,ll)
141          qfi(ij,l,:)=q(ij,ll,:)
142          Tfi(ij,l)=Temp(ij,ll)/(1+0.608*qfi(ij,l,1)) 
143       ENDDO
144    ENDDO
145   
146    IF (testcase==moist_baroclinic .OR. testcase==cyclone  ) THEN
147      CALL simple_physics(ngrid, llm, dt_phys, lat, tfi, qfi(:,:,1) , ufi, vfi, pmid, pint, pdel, 1./pdel, ps, precl, 1, .FALSE., .FALSE.)
148    ENDIF
149
150 
151    IF (testcase==moist_baroclinic .OR. testcase==cyclone .OR. testcase==supercell ) THEN
152       DO ij=1,ngrid
153          lastz=0 ;
154          DO l=1,llm
155           ll=llm+1-l
156           rho(l) = gas_constant*Temp(ij,ll) / pmid(ij,ll)
157           z(l)=lastz+ (p(ij,l)-p(ij,l+1) /g) / rho(l)
158           lastz=z(l)
159           theta(l)= Tfi(ij,l)*(1+0.608*qfi(ij,ll,1)) * pk(ij,l) / cpp
160          ENDDO
161
162          CALL KESSLER(theta(:), qfi(ij,llm:1:-1,1), qfi(ij,llm:1:-1,2), qfi(ij,llm:1:-1,3), rho(:),  &
163                       pk(ij,:), dt_phys, z(:), llm, precl(ij)) 
164         
165          DO l=1,llm
166           ll=llm+1-l
167           Tfi(ij,l) = theta(l) *(1+0.608*qfi(ij,ll,1)) * cpp / pk(ij,l)
168          ENDDO
169       ENDDO
170    ENDIF
171   
172    DO l=1,llm
173      ll=llm+1-l
174      DO ij=1,ngrid
175        CALL  tendency_terminator( lat(ij), lon(ij), qfi(ij,ll,4), qfi(ij,ll,5), dt_phys, dcl1, dcl2)
176        qfi(ij,ll,4)=qfi(ij,ll,4)+ dt_phys*dcl1
177        qfi(ij,ll,5)=qfi(ij,ll,5)+ dt_phys*dcl2
178      ENDDO
179    ENDDO
180   
181   
182    ! Mirror vertical index and compute tendencies
183    inv_dt = 1./dt_phys
184    DO l=1,llm
185       ll=llm+1-l
186       DO ij=1,ngrid
187          dTemp(ij,l) = inv_dt * ( Tfi(ij,ll)*(1+0.608*qfi(ij,ll,1)) - Temp(ij,l) )
188          du(ij,l) = inv_dt * (ufi(ij,ll) - u(ij,l))
189          dv(ij,l) = inv_dt * (vfi(ij,ll) - v(ij,l))
190          dq(ij,l,:)  = inv_dt * (qfi(ij,ll,:) - q(ij,l,:))
191       ENDDO
192    ENDDO
193
194  END SUBROUTINE compute_physics
195   
196END MODULE physics_dcmip2016_mod
197
198
Note: See TracBrowser for help on using the repository browser.