1 |
|
|
MODULE etat0dyn |
2 |
|
|
! |
3 |
|
|
!******************************************************************************* |
4 |
|
|
! Purpose: Create dynamical initial state using atmospheric fields from a |
5 |
|
|
! database of atmospheric to initialize the model. |
6 |
|
|
!------------------------------------------------------------------------------- |
7 |
|
|
! Comments: |
8 |
|
|
! |
9 |
|
|
! * This module is designed to work for Earth (and with ioipsl) |
10 |
|
|
! |
11 |
|
|
! * etat0dyn_netcdf routine can access to NetCDF data through the following |
12 |
|
|
! routine (to be called after restget): |
13 |
|
|
! CALL startget_dyn3d(varname, lon_in, lat_in, pls, workvar,& |
14 |
|
|
! champ, lon_in2, lat_in2) |
15 |
|
|
! |
16 |
|
|
! * Variables should have the following names in the NetCDF files: |
17 |
|
|
! 'U' : East ward wind (in "ECDYN.nc") |
18 |
|
|
! 'V' : Northward wind (in "ECDYN.nc") |
19 |
|
|
! 'TEMP' : Temperature (in "ECDYN.nc") |
20 |
|
|
! 'R' : Relative humidity (in "ECDYN.nc") |
21 |
|
|
! 'RELIEF' : High resolution orography (in "Relief.nc") |
22 |
|
|
! |
23 |
|
|
! * The land mask and corresponding weights can be: |
24 |
|
|
! 1) already known (in particular if etat0dyn has been called before) ; |
25 |
|
|
! in this case, ANY(masque(:,:)/=-99999.) = .TRUE. |
26 |
|
|
! 2) computed using the ocean mask from the ocean model (to ensure ocean |
27 |
|
|
! fractions are the same for atmosphere and ocean) for coupled runs. |
28 |
|
|
! File name: "o2a.nc" ; variable name: "OceMask" |
29 |
|
|
! 3) computed from topography file "Relief.nc" for forced runs. |
30 |
|
|
! |
31 |
|
|
! * There is a big mess with the longitude size. Should it be iml or iml+1 ? |
32 |
|
|
! I have chosen to use the iml+1 as an argument to this routine and we declare |
33 |
|
|
! internaly smaller fields when needed. This needs to be cleared once and for |
34 |
|
|
! all in LMDZ. A convention is required. |
35 |
|
|
!------------------------------------------------------------------------------- |
36 |
|
|
USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, histclo |
37 |
|
|
USE assert_eq_m, ONLY: assert_eq |
38 |
|
|
USE comconst_mod, ONLY: pi, cpp, kappa |
39 |
|
|
USE comvert_mod, ONLY: ap, bp, preff, pressure_exner |
40 |
|
|
USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn, itau_phy, start_time |
41 |
|
|
USE strings_mod, ONLY: strLower |
42 |
|
|
|
43 |
|
|
IMPLICIT NONE |
44 |
|
|
|
45 |
|
|
PRIVATE |
46 |
|
|
PUBLIC :: etat0dyn_netcdf |
47 |
|
|
|
48 |
|
|
include "iniprint.h" |
49 |
|
|
include "dimensions.h" |
50 |
|
|
include "paramet.h" |
51 |
|
|
include "comgeom2.h" |
52 |
|
|
include "comdissnew.h" |
53 |
|
|
REAL, SAVE :: deg2rad |
54 |
|
|
INTEGER, SAVE :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn |
55 |
|
|
REAL, ALLOCATABLE, SAVE :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:) |
56 |
|
|
CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc' |
57 |
|
|
|
58 |
|
|
CONTAINS |
59 |
|
|
|
60 |
|
|
!------------------------------------------------------------------------------- |
61 |
|
|
! |
62 |
|
|
SUBROUTINE etat0dyn_netcdf(masque, phis) |
63 |
|
|
! |
64 |
|
|
!------------------------------------------------------------------------------- |
65 |
|
|
! Purpose: Create dynamical initial states. |
66 |
|
|
!------------------------------------------------------------------------------- |
67 |
|
|
! Notes: 1) This routine is designed to work for Earth |
68 |
|
|
! 2) If masque(:,:)/=-99999., masque and phis are already known. |
69 |
|
|
! Otherwise: compute it. |
70 |
|
|
!------------------------------------------------------------------------------- |
71 |
|
|
USE control_mod |
72 |
|
|
USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz |
73 |
|
|
USE regr_pr_o3_m, ONLY: regr_pr_o3 |
74 |
|
|
USE press_coefoz_m, ONLY: press_coefoz |
75 |
|
|
USE exner_hyb_m, ONLY: exner_hyb |
76 |
|
|
USE exner_milieu_m, ONLY: exner_milieu |
77 |
|
|
USE infotrac, ONLY: nqtot, tracers |
78 |
|
|
USE filtreg_mod |
79 |
|
|
IMPLICIT NONE |
80 |
|
|
!------------------------------------------------------------------------------- |
81 |
|
|
! Arguments: |
82 |
|
|
REAL, INTENT(INOUT) :: masque(iip1,jjp1) !--- Land-ocean mask |
83 |
|
|
REAL, INTENT(INOUT) :: phis (iip1,jjp1) !--- Ground geopotential |
84 |
|
|
!------------------------------------------------------------------------------- |
85 |
|
|
! Local variables: |
86 |
|
|
CHARACTER(LEN=256) :: modname, fmt |
87 |
|
|
INTEGER :: i, j, l, ji, itau, iday, iq |
88 |
|
|
REAL :: xpn, xps, time, phystep |
89 |
|
|
REAL, DIMENSION(iip1,jjp1) :: psol |
90 |
|
|
REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d |
91 |
|
|
REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd |
92 |
|
|
REAL, DIMENSION(iip1,jjp1,llm) :: pk, pls, y, masse |
93 |
|
|
REAL, DIMENSION(iip1,jjm ,llm) :: vvent |
94 |
|
|
REAL, DIMENSION(ip1jm ,llm) :: pbarv |
95 |
|
|
REAL, DIMENSION(ip1jmp1 ,llm) :: pbaru, phi, w |
96 |
|
|
REAL, DIMENSION(ip1jmp1) :: pks |
97 |
|
|
REAL, DIMENSION(iim) :: xppn, xpps |
98 |
|
|
REAL, ALLOCATABLE :: q3d(:,:,:,:) |
99 |
|
|
!------------------------------------------------------------------------------- |
100 |
|
|
modname='etat0dyn_netcdf' |
101 |
|
|
|
102 |
|
|
deg2rad = pi/180.0 |
103 |
|
|
y(:,:,:)=0 !ym warning unitialized variable |
104 |
|
|
|
105 |
|
|
! Compute psol AND tsol, knowing phis. |
106 |
|
|
!******************************************************************************* |
107 |
|
|
CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol) |
108 |
|
|
|
109 |
|
|
! Mid-levels pressure computation |
110 |
|
|
!******************************************************************************* |
111 |
|
|
CALL pression(ip1jmp1, ap, bp, psol, p3d) !--- Update p3d |
112 |
|
|
IF(pressure_exner) THEN !--- Update pk, pks |
113 |
|
|
CALL exner_hyb (ip1jmp1,psol,p3d,pks,pk) |
114 |
|
|
ELSE |
115 |
|
|
CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk) |
116 |
|
|
END IF |
117 |
|
|
pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa) !--- Update pls |
118 |
|
|
|
119 |
|
|
! Update uvent, vvent, t3d and tpot |
120 |
|
|
!******************************************************************************* |
121 |
|
|
uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0 |
122 |
|
|
CALL startget_dyn3d('u' ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv) |
123 |
|
|
CALL startget_dyn3d('v' ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent, & |
124 |
|
|
& rlonu,rlatu(:jjm)) |
125 |
|
|
CALL startget_dyn3d('t' ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv) |
126 |
|
|
tpot(:,:,:)=t3d(:,:,:) |
127 |
|
|
CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv) |
128 |
|
|
|
129 |
|
|
WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:)) |
130 |
|
|
WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:)) |
131 |
|
|
|
132 |
|
|
! Humidity at saturation computation |
133 |
|
|
!******************************************************************************* |
134 |
|
|
WRITE(lunout,*) 'avant q_sat' |
135 |
|
|
CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat) |
136 |
|
|
WRITE(lunout,*) 'apres q_sat' |
137 |
|
|
WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:)) |
138 |
|
|
! WRITE(lunout,*) 'QSAT :',qsat(10,20,:) |
139 |
|
|
qd (:,:,:) = 0.0 |
140 |
|
|
CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv) |
141 |
|
|
ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:) |
142 |
|
|
CALL flinclo(fid_dyn) |
143 |
|
|
|
144 |
|
|
! Parameterization of ozone chemistry: |
145 |
|
|
!******************************************************************************* |
146 |
|
|
! Look for ozone tracer: |
147 |
|
|
#ifndef INCA |
148 |
|
|
DO iq=1,nqtot; IF(strLower(tracers(iq)%name)=="o3") EXIT; END DO |
149 |
|
|
IF(iq/=nqtot+1) THEN |
150 |
|
|
CALL regr_lat_time_coefoz |
151 |
|
|
CALL press_coefoz |
152 |
|
|
CALL regr_pr_o3(p3d, q3d(:,:,:,iq)) |
153 |
|
|
q3d(:,:,:,iq)=q3d(:,:,:,iq)*48./ 29. !--- Mole->mass fraction |
154 |
|
|
END IF |
155 |
|
|
#endif |
156 |
|
|
q3d(iip1,:,:,:)=q3d(1,:,:,:) |
157 |
|
|
|
158 |
|
|
! Writing |
159 |
|
|
!******************************************************************************* |
160 |
|
|
CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, & |
161 |
|
|
tetatemp, vert_prof_dissip) |
162 |
|
|
WRITE(lunout,*)'sortie inidissip' |
163 |
|
|
itau=0 |
164 |
|
|
itau_dyn=0 |
165 |
|
|
itau_phy=0 |
166 |
|
|
iday=dayref+itau/day_step |
167 |
|
|
time=FLOAT(itau-(iday-dayref)*day_step)/day_step |
168 |
|
|
IF(time>1.) THEN |
169 |
|
|
time=time-1 |
170 |
|
|
iday=iday+1 |
171 |
|
|
END IF |
172 |
|
|
day_ref=dayref |
173 |
|
|
annee_ref=anneeref |
174 |
|
|
CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi ) |
175 |
|
|
WRITE(lunout,*)'sortie geopot' |
176 |
|
|
CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, & |
177 |
|
|
phi, w, pbaru, pbarv, time+iday-dayref) |
178 |
|
|
WRITE(lunout,*)'sortie caldyn0' |
179 |
|
|
start_time = 0. |
180 |
|
|
#ifdef CPP_PARA |
181 |
|
|
CALL dynredem0_loc( "start.nc", dayref, phis) |
182 |
|
|
#else |
183 |
|
|
CALL dynredem0( "start.nc", dayref, phis) |
184 |
|
|
#endif |
185 |
|
|
WRITE(lunout,*)'sortie dynredem0' |
186 |
|
|
#ifdef CPP_PARA |
187 |
|
|
CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) |
188 |
|
|
#else |
189 |
|
|
CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) |
190 |
|
|
#endif |
191 |
|
|
WRITE(lunout,*)'sortie dynredem1' |
192 |
|
|
CALL histclo() |
193 |
|
|
|
194 |
|
|
END SUBROUTINE etat0dyn_netcdf |
195 |
|
|
! |
196 |
|
|
!------------------------------------------------------------------------------- |
197 |
|
|
|
198 |
|
|
|
199 |
|
|
!------------------------------------------------------------------------------- |
200 |
|
|
! |
201 |
|
|
SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar,& |
202 |
|
|
champ, lon_in2, lat_in2) |
203 |
|
|
!------------------------------------------------------------------------------- |
204 |
|
|
IMPLICIT NONE |
205 |
|
|
!=============================================================================== |
206 |
|
|
! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R |
207 |
|
|
! (3D fields) of file dynfname. |
208 |
|
|
!------------------------------------------------------------------------------- |
209 |
|
|
! Note: An input auxilliary field "workvar" has to be specified in two cases: |
210 |
|
|
! * for "q": the saturated humidity. |
211 |
|
|
! * for "tpot": the Exner function. |
212 |
|
|
!=============================================================================== |
213 |
|
|
! Arguments: |
214 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: var |
215 |
|
|
REAL, INTENT(IN) :: lon_in(:) ! dim (iml) |
216 |
|
|
REAL, INTENT(IN) :: lat_in(:) ! dim (jml) |
217 |
|
|
REAL, INTENT(IN) :: pls (:, :, :) ! dim (iml, jml, lml) |
218 |
|
|
REAL, INTENT(IN) :: workvar(:, :, :) ! dim (iml, jml, lml) |
219 |
|
|
REAL, INTENT(INOUT) :: champ (:, :, :) ! dim (iml, jml, lml) |
220 |
|
|
REAL, INTENT(IN) :: lon_in2(:) ! dim (iml) |
221 |
|
|
REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2) |
222 |
|
|
!------------------------------------------------------------------------------- |
223 |
|
|
! Local variables: |
224 |
|
|
CHARACTER(LEN=10) :: vname |
225 |
|
|
CHARACTER(LEN=256) :: msg, modname="startget_dyn3d" |
226 |
|
|
INTEGER :: iml, jml, jml2, lml, il |
227 |
|
|
REAL :: xppn, xpps |
228 |
|
|
!------------------------------------------------------------------------------- |
229 |
|
|
iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1), & |
230 |
|
|
& SIZE(lon_in2)], TRIM(modname)//" iml") |
231 |
|
|
jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2), & |
232 |
|
|
& TRIM(modname)//" jml") |
233 |
|
|
lml=assert_eq( SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3), & |
234 |
|
|
& TRIM(modname)//" lml") |
235 |
|
|
jml2=SIZE(lat_in2) |
236 |
|
|
|
237 |
|
|
!--- CHECK IF THE FIELD IS KNOWN |
238 |
|
|
SELECT CASE(var) |
239 |
|
|
CASE('u'); vname='U' |
240 |
|
|
CASE('v'); vname='V' |
241 |
|
|
CASE('t'); vname='TEMP' |
242 |
|
|
CASE('q'); vname='R'; msg='humidity as the saturated humidity' |
243 |
|
|
CASE('tpot'); msg='potential temperature as the Exner function' |
244 |
|
|
CASE DEFAULT; msg='No rule to extract variable '//TRIM(var) |
245 |
|
|
CALL abort_gcm(modname,TRIM(msg)//' from any data set',1) |
246 |
|
|
END SELECT |
247 |
|
|
|
248 |
|
|
!--- CHECK IF SOMETHING IS MISSING |
249 |
|
|
IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN |
250 |
|
|
msg='Could not compute '//TRIM(msg)//' is missing or constant.' |
251 |
|
|
CALL abort_gcm(modname,TRIM(msg),1) |
252 |
|
|
END IF |
253 |
|
|
|
254 |
|
|
!--- INTERPOLATE 3D FIELD IF NEEDED |
255 |
|
|
IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2, & |
256 |
|
|
lat_in2,pls,champ) |
257 |
|
|
|
258 |
|
|
!--- COMPUTE THE REQUIRED FILED |
259 |
|
|
SELECT CASE(var) |
260 |
|
|
CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO |
261 |
|
|
champ(iml,:,:)=champ(1,:,:) !--- Eastward wind |
262 |
|
|
|
263 |
|
|
CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO |
264 |
|
|
champ(iml,:,:)=champ(1,:,:) !--- Northward wind |
265 |
|
|
|
266 |
|
|
CASE('tpot','q') |
267 |
|
|
IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature |
268 |
|
|
ELSE; champ=champ*.01*workvar !--- Relative humidity |
269 |
|
|
WHERE(champ<0.) champ=1.0E-10 |
270 |
|
|
END IF |
271 |
|
|
DO il=1,lml |
272 |
|
|
xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln |
273 |
|
|
xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols |
274 |
|
|
champ(:,1 ,il) = xppn |
275 |
|
|
champ(:,jml,il) = xpps |
276 |
|
|
END DO |
277 |
|
|
END SELECT |
278 |
|
|
|
279 |
|
|
END SUBROUTINE startget_dyn3d |
280 |
|
|
! |
281 |
|
|
!------------------------------------------------------------------------------- |
282 |
|
|
|
283 |
|
|
|
284 |
|
|
!------------------------------------------------------------------------------- |
285 |
|
|
! |
286 |
|
|
SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol) |
287 |
|
|
! |
288 |
|
|
!------------------------------------------------------------------------------- |
289 |
|
|
IMPLICIT NONE |
290 |
|
|
!=============================================================================== |
291 |
|
|
! Purpose: Compute psol, knowing phis. |
292 |
|
|
!=============================================================================== |
293 |
|
|
! Arguments: |
294 |
|
|
REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) |
295 |
|
|
REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) |
296 |
|
|
REAL, INTENT(IN) :: zs (:,:) ! dim (iml,jml) |
297 |
|
|
REAL, INTENT(OUT) :: psol(:,:) ! dim (iml,jml) |
298 |
|
|
!------------------------------------------------------------------------------- |
299 |
|
|
! Local variables: |
300 |
|
|
CHARACTER(LEN=256) :: modname='start_init_dyn' |
301 |
|
|
REAL :: date, dt |
302 |
|
|
INTEGER :: iml, jml, jml2, itau(1) |
303 |
|
|
REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:,:) |
304 |
|
|
REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:) |
305 |
|
|
REAL, ALLOCATABLE :: z(:,:), ps(:,:), ts(:,:) |
306 |
|
|
!------------------------------------------------------------------------------- |
307 |
|
|
iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2), & |
308 |
|
|
& TRIM(modname)//" iml") |
309 |
|
|
jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml") |
310 |
|
|
jml2=SIZE(lat_in2) |
311 |
|
|
|
312 |
|
|
WRITE(lunout,*) 'Opening the surface analysis' |
313 |
|
|
CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn) |
314 |
|
|
WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn |
315 |
|
|
|
316 |
|
|
ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn)) |
317 |
|
|
ALLOCATE(levdyn_ini(llm_dyn)) |
318 |
|
|
CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, & |
319 |
|
|
lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn) |
320 |
|
|
|
321 |
|
|
!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS |
322 |
|
|
ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) |
323 |
|
|
lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad |
324 |
|
|
lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad |
325 |
|
|
|
326 |
|
|
ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn)) |
327 |
|
|
CALL get_var_dyn('Z',z) !--- SURFACE GEOPOTENTIAL |
328 |
|
|
CALL get_var_dyn('SP',ps) !--- SURFACE PRESSURE |
329 |
|
|
CALL get_var_dyn('ST',ts) !--- SURFACE TEMPERATURE |
330 |
|
|
! CALL flinclo(fid_dyn) |
331 |
|
|
DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) |
332 |
|
|
|
333 |
|
|
!--- PSOL IS COMPUTED IN PASCALS |
334 |
|
|
psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0 & |
335 |
|
|
& /ts(:iml-1,:)) |
336 |
|
|
psol(iml,:)=psol(1,:) |
337 |
|
|
DEALLOCATE(z,ps,ts) |
338 |
|
|
psol(:,1 )=SUM(aire(1:iml-1,1 )*psol(1:iml-1,1 ))/apoln !--- NORTH POLE |
339 |
|
|
psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols !--- SOUTH POLE |
340 |
|
|
|
341 |
|
|
CONTAINS |
342 |
|
|
|
343 |
|
|
!------------------------------------------------------------------------------- |
344 |
|
|
! |
345 |
|
|
SUBROUTINE get_var_dyn(title,field) |
346 |
|
|
! |
347 |
|
|
!------------------------------------------------------------------------------- |
348 |
|
|
USE conf_dat_m, ONLY: conf_dat2d |
349 |
|
|
IMPLICIT NONE |
350 |
|
|
!------------------------------------------------------------------------------- |
351 |
|
|
! Arguments: |
352 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: title |
353 |
|
|
REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:) |
354 |
|
|
!------------------------------------------------------------------------------- |
355 |
|
|
! Local variables: |
356 |
|
|
CHARACTER(LEN=256) :: msg |
357 |
|
|
INTEGER :: tllm |
358 |
|
|
!------------------------------------------------------------------------------- |
359 |
|
|
SELECT CASE(title) |
360 |
|
|
CASE('Z'); tllm=0; msg='geopotential' |
361 |
|
|
CASE('SP'); tllm=0; msg='surface pressure' |
362 |
|
|
CASE('ST'); tllm=llm_dyn; msg='temperature' |
363 |
|
|
END SELECT |
364 |
|
|
IF(.NOT.ALLOCATED(field)) THEN |
365 |
|
|
ALLOCATE(field(iml,jml)) |
366 |
|
|
CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana) |
367 |
|
|
CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.) |
368 |
|
|
CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana, & |
369 |
|
|
lon_in, lat_in, lon_in2, lat_in2, field) |
370 |
|
|
ELSE IF(SIZE(field)/=SIZE(z)) THEN |
371 |
|
|
msg='The '//TRIM(msg)//' field we have does not have the right size' |
372 |
|
|
CALL abort_gcm(TRIM(modname),msg,1) |
373 |
|
|
END IF |
374 |
|
|
|
375 |
|
|
END SUBROUTINE get_var_dyn |
376 |
|
|
! |
377 |
|
|
!------------------------------------------------------------------------------- |
378 |
|
|
|
379 |
|
|
END SUBROUTINE start_init_dyn |
380 |
|
|
! |
381 |
|
|
!------------------------------------------------------------------------------- |
382 |
|
|
|
383 |
|
|
|
384 |
|
|
!------------------------------------------------------------------------------- |
385 |
|
|
! |
386 |
|
|
SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d) |
387 |
|
|
! |
388 |
|
|
!------------------------------------------------------------------------------- |
389 |
|
|
USE conf_dat_m, ONLY: conf_dat3d |
390 |
|
|
USE pchsp_95_m, ONLY: pchsp_95 |
391 |
|
|
USE pchfe_95_m, ONLY: pchfe_95 |
392 |
|
|
IMPLICIT NONE |
393 |
|
|
!------------------------------------------------------------------------------- |
394 |
|
|
! Arguments: |
395 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: var |
396 |
|
|
REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) |
397 |
|
|
REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) |
398 |
|
|
REAL, INTENT(IN) :: pls_in(:,:,:) ! dim (iml,jml,lml) |
399 |
|
|
REAL, INTENT(OUT) :: var3d (:,:,:) ! dim (iml,jml,lml) |
400 |
|
|
!------------------------------------------------------------------------------- |
401 |
|
|
! Local variables: |
402 |
|
|
CHARACTER(LEN=256) :: modname='start_inter_3d' |
403 |
|
|
LOGICAL :: skip |
404 |
|
|
REAL :: chmin, chmax |
405 |
|
|
INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr |
406 |
|
|
INTEGER :: n_extrap ! Extrapolated points number |
407 |
|
|
REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:) |
408 |
|
|
REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:) |
409 |
|
|
REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:) |
410 |
|
|
!------------------------------------------------------------------------------- |
411 |
|
|
iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml") |
412 |
|
|
jml=assert_eq(SIZE(lat_in), SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml") |
413 |
|
|
lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2) |
414 |
|
|
|
415 |
|
|
WRITE(lunout, *)'Going into flinget to extract the 3D field.' |
416 |
|
|
IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn)) |
417 |
|
|
CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d) |
418 |
|
|
|
419 |
|
|
!--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS |
420 |
|
|
ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn)) |
421 |
|
|
lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad |
422 |
|
|
lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad |
423 |
|
|
|
424 |
|
|
!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS |
425 |
|
|
ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn)) |
426 |
|
|
CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, & |
427 |
|
|
lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.) |
428 |
|
|
DEALLOCATE(lon_ini, lat_ini) |
429 |
|
|
|
430 |
|
|
!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro |
431 |
|
|
ALLOCATE(var_tmp3d(iml,jml,llm_dyn)) |
432 |
|
|
DO il = 1,llm_dyn |
433 |
|
|
CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il), & |
434 |
|
|
lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il)) |
435 |
|
|
END DO |
436 |
|
|
DEALLOCATE(lon_rad, lat_rad) |
437 |
|
|
|
438 |
|
|
!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND |
439 |
|
|
ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn)) |
440 |
|
|
ax = lev_dyn(llm_dyn:1:-1) |
441 |
|
|
skip = .FALSE. |
442 |
|
|
n_extrap = 0 |
443 |
|
|
DO ij=1, jml |
444 |
|
|
DO ii=1, iml-1 |
445 |
|
|
ay = var_tmp3d(ii, ij, llm_dyn:1:-1) |
446 |
|
|
yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.) |
447 |
|
|
CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), & |
448 |
|
|
var3d(ii, ij, lml:1:-1), ierr) |
449 |
|
|
IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1) |
450 |
|
|
n_extrap = n_extrap + ierr |
451 |
|
|
END DO |
452 |
|
|
END DO |
453 |
|
|
IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap |
454 |
|
|
var3d(iml, :, :) = var3d(1, :, :) |
455 |
|
|
|
456 |
|
|
DO il=1, lml |
457 |
|
|
CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax) |
458 |
|
|
WRITE(lunout, *)' '//TRIM(var)//' min max l ', il, chmin, chmax |
459 |
|
|
END DO |
460 |
|
|
|
461 |
|
|
END SUBROUTINE start_inter_3d |
462 |
|
|
! |
463 |
|
|
!------------------------------------------------------------------------------- |
464 |
|
|
|
465 |
|
|
|
466 |
|
|
!------------------------------------------------------------------------------- |
467 |
|
|
! |
468 |
|
|
SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo) |
469 |
|
|
! |
470 |
|
|
!------------------------------------------------------------------------------- |
471 |
|
|
USE inter_barxy_m, ONLY: inter_barxy |
472 |
|
|
IMPLICIT NONE |
473 |
|
|
!------------------------------------------------------------------------------- |
474 |
|
|
! Arguments: |
475 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: nam |
476 |
|
|
LOGICAL, INTENT(IN) :: ibeg |
477 |
|
|
REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) |
478 |
|
|
REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) |
479 |
|
|
REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1) |
480 |
|
|
REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2) |
481 |
|
|
REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1) |
482 |
|
|
!------------------------------------------------------------------------------- |
483 |
|
|
! Local variables: |
484 |
|
|
CHARACTER(LEN=256) :: modname="interp_startvar" |
485 |
|
|
INTEGER :: ii, jj, i1, j1, j2 |
486 |
|
|
REAL, ALLOCATABLE :: vtmp(:,:) |
487 |
|
|
!------------------------------------------------------------------------------- |
488 |
|
|
ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii") |
489 |
|
|
jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj") |
490 |
|
|
i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1") |
491 |
|
|
j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1") |
492 |
|
|
j2=SIZE(lat2) |
493 |
|
|
ALLOCATE(vtmp(i1-1,j1)) |
494 |
|
|
IF(ibeg.AND.prt_level>1) THEN |
495 |
|
|
WRITE(lunout,*)"---------------------------------------------------------" |
496 |
|
|
WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" |
497 |
|
|
WRITE(lunout,*)"---------------------------------------------------------" |
498 |
|
|
END IF |
499 |
|
|
CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) |
500 |
|
|
CALL gr_int_dyn(vtmp, varo, i1-1, j1) |
501 |
|
|
|
502 |
|
|
END SUBROUTINE interp_startvar |
503 |
|
|
! |
504 |
|
|
!------------------------------------------------------------------------------- |
505 |
|
|
|
506 |
|
|
END MODULE etat0dyn |
507 |
|
|
! |
508 |
|
|
!******************************************************************************* |