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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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