/[lmdze]/trunk/libf/dyn3d/etat0.f90
ViewVC logotype

Annotation of /trunk/libf/dyn3d/etat0.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 36 - (hide annotations)
Thu Dec 2 17:11:04 2010 UTC (13 years, 5 months ago) by guez
File size: 11829 byte(s)
Now using the library "NR_util".

1 guez 3 module etat0_mod
2    
3     use indicesol, only: nbsrf
4     use dimphy, only: klon
5    
6     IMPLICIT NONE
7    
8     REAL pctsrf(klon, nbsrf)
9 guez 15 ! ("pctsrf(i, :)" is the composition of the surface at horizontal
10     ! position "i")
11 guez 3
12     private nbsrf, klon
13    
14     contains
15    
16     SUBROUTINE etat0
17    
18     ! From "etat0_netcdf.F", version 1.3 2005/05/25 13:10:09
19    
20 guez 15 ! This subroutine creates "mask".
21 guez 3
22 guez 27 use caldyn0_m, only: caldyn0
23 guez 3 use comconst, only: dtvr, daysec, cpp, kappa, pi
24 guez 27 use comgeom, only: rlatu, rlonv, rlonu, rlatv, aire_2d, apoln, apols, &
25     cu_2d, cv_2d
26 guez 3 use comvert, only: ap, bp, preff, pa
27 guez 27 use conf_gcm_m, only: day_step, iphysiq, dayref, anneeref
28     use dimens_m, only: iim, jjm, llm, nqmx
29 guez 3 use dimphy, only: zmasq
30     use dimsoil, only: nsoilmx
31 guez 27 use dynredem0_m, only: dynredem0
32     use dynredem1_m, only: dynredem1
33     use exner_hyb_m, only: exner_hyb
34 guez 3 use grid_atob, only: grille_m
35     use grid_change, only: init_dyn_phy, dyn_phy
36 guez 27 use indicesol, only: is_oce, is_sic, is_ter, is_lic, epsfra
37 guez 18 use iniadvtrac_m, only: iniadvtrac
38 guez 27 use inidissip_m, only: inidissip
39     use inigeom_m, only: inigeom
40 guez 32 USE flincom, only: flinclo, flinopen_nozoom, flininfo
41     use flinget_m, only: flinget
42 guez 30 use histcom, only: histclo
43 guez 27 use paramet_m, only: ip1jm, ip1jmp1
44     use phyredem_m, only: phyredem
45 guez 10 use pressure_var, only: pls, p3d
46 guez 27 use q_sat_m, only: q_sat
47 guez 7 use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
48     use regr_pr_o3_m, only: regr_pr_o3
49 guez 27 use serre, only: alphax
50     USE start_init_orog_m, only: start_init_orog, mask, phis
51     use start_init_phys_m, only: qsol_2d
52     use startdyn, only: start_inter_3d, start_init_dyn
53 guez 28 use temps, only: itau_phy, annee_ref, day_ref
54 guez 3
55     ! Variables local to the procedure:
56    
57     REAL latfi(klon), lonfi(klon)
58     ! (latitude and longitude of a point of the scalar grid identified
59     ! by a simple index, in °)
60    
61     REAL, dimension(iim + 1, jjm + 1, llm):: uvent, t3d, tpot
62     REAL vvent(iim + 1, jjm, llm)
63    
64     REAL q3d(iim + 1, jjm + 1, llm, nqmx)
65     ! (mass fractions of trace species
66     ! "q3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
67     ! and pressure level "pls(i, j, l)".)
68    
69     real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
70     REAL tsol(klon), qsol(klon), sn(klon)
71     REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf)
72     REAL albe(klon, nbsrf), evap(klon, nbsrf)
73     REAL alblw(klon, nbsrf)
74     REAL tsoil(klon, nsoilmx, nbsrf)
75     REAL radsol(klon), rain_fall(klon), snow_fall(klon)
76     REAL solsw(klon), sollw(klon), fder(klon)
77     !IM "slab" ocean
78     REAL tslab(klon)
79     real seaice(klon) ! kg m-2
80     REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
81     REAL rugmer(klon)
82     real, dimension(iim + 1, jjm + 1):: relief, zstd_2d, zsig_2d, zgam_2d
83     real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
84     real, dimension(iim + 1, jjm + 1):: tsol_2d, psol
85     REAL zmea(klon), zstd(klon)
86     REAL zsig(klon), zgam(klon)
87     REAL zthe(klon)
88     REAL zpic(klon), zval(klon)
89     REAL t_ancien(klon, llm), q_ancien(klon, llm) !
90     REAL run_off_lic_0(klon)
91     real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
92     ! déclarations pour lecture glace de mer
93     INTEGER iml_lic, jml_lic, llm_tmp, ttm_tmp
94     INTEGER itaul(1), fid
95     REAL lev(1), date
96     REAL, ALLOCATABLE:: lon_lic(:, :), lat_lic(:, :)
97     REAL, ALLOCATABLE:: dlon_lic(:), dlat_lic(:)
98     REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
99     REAL flic_tmp(iim + 1, jjm + 1) !fraction land ice temporary
100    
101     INTEGER l, ji
102    
103     REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
104     real pks(iim + 1, jjm + 1)
105    
106     REAL masse(iim + 1, jjm + 1, llm)
107     REAL phi(ip1jmp1, llm)
108     REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
109     REAL w(ip1jmp1, llm)
110     REAL phystep
111 guez 28 real trash
112 guez 3
113     !---------------------------------
114    
115     print *, "Call sequence information: etat0"
116    
117     dtvr = daysec / real(day_step)
118     print *, 'dtvr = ', dtvr
119    
120 guez 36 ! Construct a grid:
121    
122 guez 3 pa = 5e4
123     CALL iniconst
124     CALL inigeom
125     CALL inifilr
126    
127     latfi(1) = 90.
128     latfi(2:klon-1) = pack(spread(rlatu(2:jjm), 1, iim), .true.) * 180. / pi
129     ! (with conversion to degrees)
130     latfi(klon) = - 90.
131    
132     lonfi(1) = 0.
133     lonfi(2:klon-1) = pack(spread(rlonv(:iim), 2, jjm - 1), .true.) * 180. / pi
134     ! (with conversion to degrees)
135     lonfi(klon) = 0.
136    
137     call start_init_orog(relief, zstd_2d, zsig_2d, zgam_2d, zthe_2d, zpic_2d, &
138 guez 15 zval_2d) ! also compute "mask" and "phis"
139 guez 3 call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
140 guez 15 zmasq = pack(mask, dyn_phy)
141 guez 3 PRINT *, 'Masque construit'
142    
143     CALL start_init_dyn(tsol_2d, psol) ! also compute "qsol_2d"
144    
145     ! Compute pressure on intermediate levels:
146 guez 15 forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol
147 guez 3 CALL exner_hyb(psol, p3d, pks, pk)
148     IF (MINVAL(pk) == MAXVAL(pk)) stop '"pk" should not be constant'
149    
150 guez 36 pls = preff * (pk / cpp)**(1. / kappa)
151     PRINT *, "minval(pls) = ", minval(pls)
152     print *, "maxval(pls) = ", maxval(pls)
153 guez 3
154 guez 23 call start_inter_3d('U', rlonv, rlatv, pls, uvent)
155 guez 3 forall (l = 1: llm) uvent(:iim, :, l) = uvent(:iim, :, l) * cu_2d(:iim, :)
156     uvent(iim+1, :, :) = uvent(1, :, :)
157    
158 guez 23 call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vvent)
159 guez 3 forall (l = 1: llm) vvent(:iim, :, l) = vvent(:iim, :, l) * cv_2d(:iim, :)
160     vvent(iim + 1, :, :) = vvent(1, :, :)
161    
162 guez 23 call start_inter_3d('TEMP', rlonu, rlatv, pls, t3d)
163 guez 36 PRINT *, 'minval(t3d) = ', minval(t3d)
164     print *, "maxval(t3d) = ", maxval(t3d)
165 guez 3
166     tpot(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
167     tpot(iim + 1, :, :) = tpot(1, :, :)
168     DO l=1, llm
169     tpot(:, 1, l) = SUM(aire_2d(:, 1) * tpot(:, 1, l)) / apoln
170     tpot(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * tpot(:, jjm + 1, l)) &
171     / apols
172     ENDDO
173    
174     ! Calcul de l'humidité à saturation :
175 guez 36 qsat = q_sat(t3d, pls)
176     PRINT *, "minval(qsat) = ", minval(qsat)
177     print *, "maxval(qsat) = ", maxval(qsat)
178 guez 3 IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'
179    
180     ! Water vapor:
181 guez 23 call start_inter_3d('R', rlonu, rlatv, pls, q3d(:, :, :, 1))
182     q3d(:, :, :, 1) = 0.01 * q3d(:, :, :, 1) * qsat
183 guez 3 WHERE (q3d(:, :, :, 1) < 0.) q3d(:, :, :, 1) = 1E-10
184     DO l = 1, llm
185     q3d(:, 1, l, 1) = SUM(aire_2d(:, 1) * q3d(:, 1, l, 1)) / apoln
186     q3d(:, jjm + 1, l, 1) &
187     = SUM(aire_2d(:, jjm + 1) * q3d(:, jjm + 1, l, 1)) / apols
188     ENDDO
189    
190     q3d(:, :, :, 2:4) = 0. ! liquid water, radon and lead
191    
192 guez 5 if (nqmx >= 5) then
193     ! Ozone:
194 guez 7 call regr_lat_time_coefoz
195     call regr_pr_o3(q3d(:, :, :, 5))
196     ! Convert from mole fraction to mass fraction:
197     q3d(:, :, :, 5) = q3d(:, :, :, 5) * 48. / 29.
198 guez 5 end if
199 guez 3
200 guez 13 tsol = pack(tsol_2d, dyn_phy)
201     qsol = pack(qsol_2d, dyn_phy)
202     sn = 0. ! snow
203     radsol = 0.
204     tslab = 0. ! IM "slab" ocean
205     seaice = 0.
206     rugmer = 0.001
207     zmea = pack(relief, dyn_phy)
208     zstd = pack(zstd_2d, dyn_phy)
209     zsig = pack(zsig_2d, dyn_phy)
210     zgam = pack(zgam_2d, dyn_phy)
211     zthe = pack(zthe_2d, dyn_phy)
212     zpic = pack(zpic_2d, dyn_phy)
213     zval = pack(zval_2d, dyn_phy)
214 guez 3
215     ! On initialise les sous-surfaces:
216     ! Lecture du fichier glace de terre pour fixer la fraction de terre
217     ! et de glace de terre:
218     CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, &
219     ttm_tmp, fid)
220     ALLOCATE(lat_lic(iml_lic, jml_lic))
221     ALLOCATE(lon_lic(iml_lic, jml_lic))
222     ALLOCATE(dlon_lic(iml_lic))
223     ALLOCATE(dlat_lic(jml_lic))
224     ALLOCATE(fraclic(iml_lic, jml_lic))
225 guez 32 CALL flinopen_nozoom(iml_lic, jml_lic, &
226 guez 28 llm_tmp, lon_lic, lat_lic, lev, ttm_tmp, itaul, date, trash, &
227 guez 3 fid)
228     CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp &
229     , 1, 1, fraclic)
230     CALL flinclo(fid)
231    
232     ! Interpolation sur la grille T du modèle :
233     PRINT *, 'Dimensions de "landice"'
234     print *, "iml_lic = ", iml_lic
235     print *, "jml_lic = ", jml_lic
236    
237     ! Si les coordonnées sont en degrés, on les transforme :
238 guez 15 IF (MAXVAL( lon_lic ) > pi) THEN
239     lon_lic = lon_lic * pi / 180.
240 guez 3 ENDIF
241 guez 15 IF (maxval( lat_lic ) > pi) THEN
242     lat_lic = lat_lic * pi/ 180.
243 guez 3 ENDIF
244    
245 guez 13 dlon_lic = lon_lic(:, 1)
246     dlat_lic = lat_lic(1, :)
247 guez 3
248     flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
249     rlatu)
250     flic_tmp(iim + 1, :) = flic_tmp(1, :)
251    
252     ! Passage sur la grille physique
253 guez 15 pctsrf = 0.
254 guez 3 pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
255     ! Adéquation avec le maque terre/mer
256     WHERE (pctsrf(:, is_lic) < EPSFRA ) pctsrf(:, is_lic) = 0.
257 guez 13 WHERE (zmasq < EPSFRA) pctsrf(:, is_lic) = 0.
258     pctsrf(:, is_ter) = zmasq
259     where (zmasq > EPSFRA)
260     where (pctsrf(:, is_lic) >= zmasq)
261     pctsrf(:, is_lic) = zmasq
262 guez 3 pctsrf(:, is_ter) = 0.
263     elsewhere
264 guez 13 pctsrf(:, is_ter) = zmasq - pctsrf(:, is_lic)
265 guez 3 where (pctsrf(:, is_ter) < EPSFRA)
266     pctsrf(:, is_ter) = 0.
267 guez 13 pctsrf(:, is_lic) = zmasq
268 guez 3 end where
269     end where
270     end where
271    
272     ! Sous-surface océan et glace de mer (pour démarrer on met glace
273     ! de mer à 0) :
274 guez 13 pctsrf(:, is_oce) = 1. - zmasq
275 guez 3 WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
276    
277     ! Vérification que somme des sous-surfaces vaut 1:
278 guez 15 ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
279     IF (ji /= 0) then
280     PRINT *, 'Problème répartition sous maille pour ', ji, 'points'
281     end IF
282 guez 3
283     ! Calcul intermédiaire:
284     CALL massdair(p3d, masse)
285    
286     print *, 'ALPHAX = ', alphax
287    
288     forall (l = 1:llm)
289     masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
290     masse(:, jjm + 1, l) = &
291     SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
292     END forall
293    
294     ! Initialisation pour traceurs:
295 guez 5 call iniadvtrac
296 guez 27 CALL inidissip
297 guez 3 itau_phy = 0
298     day_ref = dayref
299     annee_ref = anneeref
300    
301 guez 22 CALL geopot(ip1jmp1, tpot, pk , pks, phis, phi)
302 guez 23 CALL caldyn0(uvent, vvent, tpot, psol, masse, pk, phis, phi, w, pbaru, &
303     pbarv)
304 guez 5 CALL dynredem0("start.nc", dayref, phis)
305 guez 27 CALL dynredem1("start.nc", vvent, uvent, tpot, q3d, masse, psol, itau=0)
306 guez 3
307     ! Ecriture état initial physique:
308     print *, "iphysiq = ", iphysiq
309     phystep = dtvr * REAL(iphysiq)
310     print *, 'phystep = ', phystep
311    
312     ! Initialisations :
313     tsolsrf(:, is_ter) = tsol
314     tsolsrf(:, is_lic) = tsol
315     tsolsrf(:, is_oce) = tsol
316     tsolsrf(:, is_sic) = tsol
317     snsrf(:, is_ter) = sn
318     snsrf(:, is_lic) = sn
319     snsrf(:, is_oce) = sn
320     snsrf(:, is_sic) = sn
321     albe(:, is_ter) = 0.08
322     albe(:, is_lic) = 0.6
323     albe(:, is_oce) = 0.5
324     albe(:, is_sic) = 0.6
325     alblw = albe
326 guez 15 evap = 0.
327 guez 3 qsolsrf(:, is_ter) = 150.
328     qsolsrf(:, is_lic) = 150.
329     qsolsrf(:, is_oce) = 150.
330     qsolsrf(:, is_sic) = 150.
331     tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
332     rain_fall = 0.
333     snow_fall = 0.
334     solsw = 165.
335     sollw = -53.
336     t_ancien = 273.15
337     q_ancien = 0.
338     agesno = 0.
339     !IM "slab" ocean
340 guez 13 tslab = tsolsrf(:, is_oce)
341 guez 3 seaice = 0.
342    
343 guez 13 frugs(:, is_oce) = rugmer
344     frugs(:, is_ter) = MAX(1.e-05, zstd * zsig / 2)
345     frugs(:, is_lic) = MAX(1.e-05, zstd * zsig / 2)
346 guez 3 frugs(:, is_sic) = 0.001
347     fder = 0.
348     clwcon = 0.
349     rnebcon = 0.
350     ratqs = 0.
351     run_off_lic_0 = 0.
352    
353 guez 13 call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
354 guez 3 tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
355     evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
356 guez 13 agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
357 guez 3 t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
358     CALL histclo
359    
360     END SUBROUTINE etat0
361    
362     end module etat0_mod

  ViewVC Help
Powered by ViewVC 1.1.21