/[lmdze]/trunk/Sources/dyn3d/etat0.f
ViewVC logotype

Annotation of /trunk/Sources/dyn3d/etat0.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 30 - (hide annotations)
Thu Apr 1 09:07:28 2010 UTC (14 years, 1 month ago) by guez
Original Path: trunk/libf/dyn3d/etat0.f90
File size: 11957 byte(s)
Imported Source files of the external library "IOIPSL_Lionel" into
"libf/IOIPSL".

Split "cray.f90" into "scopy.f90" and "ssum.f90".

Rewrote "leapfrog" in order to have a clearer algorithmic structure.

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

  ViewVC Help
Powered by ViewVC 1.1.21