/[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 10 by guez, Fri Apr 18 14:45:53 2008 UTC trunk/dyn3d/etat0.f revision 97 by guez, Fri Apr 25 14:58:31 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 q_sat_m, only: q_sat      use q_sat_m, only: q_sat
     use exner_hyb_m, only: exner_hyb  
     use advtrac_m, only: iniadvtrac  
     use pressure_var, only: pls, p3d  
     use dynredem0_m, only: dynredem0  
47      use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz      use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
48      use regr_pr_o3_m, only: regr_pr_o3      use regr_pr_o3_m, only: regr_pr_o3
49        use serre, only: alphax
50        use startdyn, only: start_init_dyn
51        USE start_init_orog_m, only: start_init_orog, mask
52        use start_init_phys_m, only: start_init_phys
53        use start_inter_3d_m, only: start_inter_3d
54        use temps, only: itau_phy, annee_ref, day_ref
55    
56      ! Variables local to the procedure:      ! Variables local to the procedure:
57    
58      REAL latfi(klon), lonfi(klon)      REAL latfi(klon), lonfi(klon)
59      ! (latitude and longitude of a point of the scalar grid identified      ! (latitude and longitude of a point of the scalar grid identified
60      ! by a simple index, in °)      ! by a simple index, in degrees)
61    
62      REAL, dimension(iim + 1, jjm + 1, llm):: uvent, t3d, tpot      REAL, dimension(iim + 1, jjm + 1, llm):: ucov, t3d, teta
63      REAL vvent(iim + 1, jjm, llm)      REAL vcov(iim + 1, jjm, llm)
64    
65      REAL q3d(iim + 1, jjm + 1, llm, nqmx)      REAL q(iim + 1, jjm + 1, llm, nqmx)
66      ! (mass fractions of trace species      ! (mass fractions of trace species
67      ! "q3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"      ! "q(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
68      ! and pressure level "pls(i, j, l)".)      ! and pressure level "pls(i, j, l)".)
69    
70      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 80  contains
80      real seaice(klon) ! kg m-2      real seaice(klon) ! kg m-2
81      REAL frugs(klon, nbsrf), agesno(klon, nbsrf)      REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
82      REAL rugmer(klon)      REAL rugmer(klon)
83      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
84        real, dimension(iim + 1, jjm + 1):: zmea_2d, zstd_2d, zsig_2d, zgam_2d
85      real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d      real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
86      real, dimension(iim + 1, jjm + 1):: tsol_2d, psol      real, dimension(iim + 1, jjm + 1):: tsol_2d, qsol_2d, ps
87      REAL zmea(klon), zstd(klon)      REAL zmea(klon), zstd(klon)
88      REAL zsig(klon), zgam(klon)      REAL zsig(klon), zgam(klon)
89      REAL zthe(klon)      REAL zthe(klon)
90      REAL zpic(klon), zval(klon)      REAL zpic(klon), zval(klon)
91      REAL rugsrel(klon)      REAL t_ancien(klon, llm), q_ancien(klon, llm)
     REAL t_ancien(klon, llm), q_ancien(klon, llm)      !  
92      REAL run_off_lic_0(klon)      REAL run_off_lic_0(klon)
93      real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)      real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
94      ! déclarations pour lecture glace de mer  
95      INTEGER iml_lic, jml_lic, llm_tmp, ttm_tmp      ! D\'eclarations pour lecture glace de mer :
96      INTEGER itaul(1), fid      INTEGER iml_lic, jml_lic
97      REAL lev(1), date      INTEGER ncid, varid
98      REAL, ALLOCATABLE:: lon_lic(:, :), lat_lic(:, :)      REAL, pointer:: dlon_lic(:), dlat_lic(:)
     REAL, ALLOCATABLE:: dlon_lic(:), dlat_lic(:)  
99      REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice      REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
100      REAL flic_tmp(iim + 1, jjm + 1) !fraction land ice temporary      REAL flic_tmp(iim + 1, jjm + 1) ! fraction land ice temporary
101    
102      INTEGER l, ji      INTEGER l, ji
103    
# Line 100  contains Line 105  contains
105      real pks(iim + 1, jjm + 1)      real pks(iim + 1, jjm + 1)
106    
107      REAL masse(iim + 1, jjm + 1, llm)      REAL masse(iim + 1, jjm + 1, llm)
108      REAL phi(ip1jmp1, llm)      REAL phi(iim + 1, jjm + 1, llm)
109      REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)      REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
110      REAL w(ip1jmp1, llm)      REAL w(iim + 1, jjm + 1, llm)
111      REAL phystep  
112      INTEGER radpas      real sig1(klon, llm) ! section adiabatic updraft
113        real w01(klon, llm) ! vertical velocity within adiabatic updraft
114    
115        real pls(iim + 1, jjm + 1, llm)
116        ! (pressure at mid-layer of LMDZ grid, in Pa)
117        ! "pls(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
118        ! for layer "l")
119    
120        REAL p3d(iim + 1, jjm + 1, llm+1) ! pressure at layer interfaces, in Pa
121        ! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
122        ! for interface "l")
123    
124      !---------------------------------      !---------------------------------
125    
126      print *, "Call sequence information: etat0"      print *, "Call sequence information: etat0"
127    
128      ! Construct a grid:      CALL iniconst
129    
130      dtvr = daysec / real(day_step)      ! Construct a grid:
     print *, 'dtvr = ', dtvr  
131    
132      pa = 5e4      pa = 5e4
133      CALL iniconst      CALL disvert
134      CALL inigeom      CALL inigeom
135      CALL inifilr      CALL inifilr
136    
# Line 130  contains Line 144  contains
144      ! (with conversion to degrees)      ! (with conversion to degrees)
145      lonfi(klon) = 0.      lonfi(klon) = 0.
146    
147      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, &
148           zval_2d) ! also compute "masque" and "phis"           zpic_2d, zval_2d) ! also compute "mask"
149      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
150      zmasq = pack(masque, dyn_phy)      zmasq = pack(mask, dyn_phy)
151      PRINT *, 'Masque construit'      PRINT *, 'Masque construit'
152    
153      CALL start_init_dyn(tsol_2d, psol) ! also compute "qsol_2d"      call start_init_phys(tsol_2d, qsol_2d)
154        CALL start_init_dyn(tsol_2d, phis, ps)
155    
156      ! Compute pressure on intermediate levels:      ! Compute pressure on intermediate levels:
157      forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol(:, :)      forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
158      CALL exner_hyb(psol, p3d, pks, pk)      CALL exner_hyb(ps, p3d, pks, pk)
159      IF (MINVAL(pk) == MAXVAL(pk)) stop '"pk" should not be constant'      call assert(MINVAL(pk) /= MAXVAL(pk), '"pk" should not be constant')
160    
161      pls(:, :, :) = preff * (pk(:, :, :) / cpp)**(1. / kappa)      pls = preff * (pk / cpp)**(1. / kappa)
162      PRINT *, "minval(pls(:, :, :)) = ", minval(pls(:, :, :))      PRINT *, "minval(pls) = ", minval(pls)
163      print *, "maxval(pls(:, :, :)) = ", maxval(pls(:, :, :))      print *, "maxval(pls) = ", maxval(pls)
164    
165      uvent(:, :, :) = start_inter_3d('U', rlonv, rlatv, pls)      call start_inter_3d('U', rlonv, rlatv, pls, ucov)
166      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, :)
167      uvent(iim+1, :, :) = uvent(1, :, :)      ucov(iim+1, :, :) = ucov(1, :, :)
168    
169      vvent(:, :, :) = start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :))      call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vcov)
170      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, :)
171      vvent(iim + 1, :, :) = vvent(1, :, :)      vcov(iim + 1, :, :) = vcov(1, :, :)
172    
173      t3d(:, :, :) = start_inter_3d('TEMP', rlonu, rlatv, pls)      call start_inter_3d('TEMP', rlonu, rlatv, pls, t3d)
174      PRINT *,  'minval(t3d(:, :, :)) = ', minval(t3d(:, :, :))      PRINT *, 'minval(t3d) = ', minval(t3d)
175      print *, "maxval(t3d(:, :, :)) = ", maxval(t3d(:, :, :))      print *, "maxval(t3d) = ", maxval(t3d)
176    
177      tpot(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)      teta(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
178      tpot(iim + 1, :, :) = tpot(1, :, :)      teta(iim + 1, :, :) = teta(1, :, :)
179      DO l=1, llm      DO l = 1, llm
180         tpot(:, 1, l) = SUM(aire_2d(:, 1) * tpot(:, 1, l)) / apoln         teta(:, 1, l) = SUM(aire_2d(:, 1) * teta(:, 1, l)) / apoln
181         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)) &
182              / apols              / apols
183      ENDDO      ENDDO
184    
185      ! Calcul de l'humidité à saturation :      ! Calcul de l'humidit\'e \`a saturation :
186      qsat(:, :, :) = q_sat(t3d, pls)      qsat = q_sat(t3d, pls)
187      PRINT *, "minval(qsat(:, :, :)) = ", minval(qsat(:, :, :))      PRINT *, "minval(qsat) = ", minval(qsat)
188      print *, "maxval(qsat(:, :, :)) = ", maxval(qsat(:, :, :))      print *, "maxval(qsat) = ", maxval(qsat)
189      IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'      IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'
190    
191      ! Water vapor:      ! Water vapor:
192      q3d(:, :, :, 1) = 0.01 * start_inter_3d('R', rlonu, rlatv, pls) * qsat      call start_inter_3d('R', rlonu, rlatv, pls, q(:, :, :, 1))
193      WHERE (q3d(:, :, :, 1) < 0.) q3d(:, :, :, 1) = 1E-10      q(:, :, :, 1) = 0.01 * q(:, :, :, 1) * qsat
194        WHERE (q(:, :, :, 1) < 0.) q(:, :, :, 1) = 1E-10
195      DO l = 1, llm      DO l = 1, llm
196         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
197         q3d(:, jjm + 1, l, 1) &         q(:, jjm + 1, l, 1) &
198              = SUM(aire_2d(:, jjm + 1) * q3d(:, jjm + 1, l, 1)) / apols              = SUM(aire_2d(:, jjm + 1) * q(:, jjm + 1, l, 1)) / apols
199      ENDDO      ENDDO
200    
201      q3d(:, :, :, 2:4) = 0. ! liquid water, radon and lead      q(:, :, :, 2:4) = 0. ! liquid water, radon and lead
202    
203      if (nqmx >= 5) then      if (nqmx >= 5) then
204         ! Ozone:         ! Ozone:
205         call regr_lat_time_coefoz         call regr_lat_time_coefoz
206         call regr_pr_o3(q3d(:, :, :, 5))         call regr_pr_o3(p3d, q(:, :, :, 5))
207         ! Convert from mole fraction to mass fraction:         ! Convert from mole fraction to mass fraction:
208         q3d(:, :, :, 5) = q3d(:, :, :, 5)  * 48. / 29.         q(:, :, :, 5) = q(:, :, :, 5) * 48. / 29.
209      end if      end if
210    
211      tsol(:) = pack(tsol_2d, dyn_phy)      tsol = pack(tsol_2d, dyn_phy)
212      qsol(:) = pack(qsol_2d, dyn_phy)      qsol = pack(qsol_2d, dyn_phy)
213      sn(:) = 0. ! snow      sn = 0. ! snow
214      radsol(:) = 0.      radsol = 0.
215      tslab(:) = 0. ! IM "slab" ocean      tslab = 0. ! IM "slab" ocean
216      seaice(:) = 0.      seaice = 0.
217      rugmer(:) = 0.001      rugmer = 0.001
218      zmea(:) = pack(relief, dyn_phy)      zmea = pack(zmea_2d, dyn_phy)
219      zstd(:) = pack(zstd_2d, dyn_phy)      zstd = pack(zstd_2d, dyn_phy)
220      zsig(:) = pack(zsig_2d, dyn_phy)      zsig = pack(zsig_2d, dyn_phy)
221      zgam(:) = pack(zgam_2d, dyn_phy)      zgam = pack(zgam_2d, dyn_phy)
222      zthe(:) = pack(zthe_2d, dyn_phy)      zthe = pack(zthe_2d, dyn_phy)
223      zpic(:) = pack(zpic_2d, dyn_phy)      zpic = pack(zpic_2d, dyn_phy)
224      zval(:) = pack(zval_2d, dyn_phy)      zval = pack(zval_2d, dyn_phy)
   
     rugsrel(:) = 0.  
     IF (ok_orodr) rugsrel(:) = MAX(1.e-05, zstd(:) * zsig(:) / 2)  
225    
226      ! On initialise les sous-surfaces:      ! On initialise les sous-surfaces.
227      ! Lecture du fichier glace de terre pour fixer la fraction de terre      ! Lecture du fichier glace de terre pour fixer la fraction de terre
228      ! et de glace de terre:      ! et de glace de terre :
229      CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, &  
230           ttm_tmp, fid)      call nf95_open("landiceref.nc", nf90_nowrite, ncid)
231      ALLOCATE(lat_lic(iml_lic, jml_lic))  
232      ALLOCATE(lon_lic(iml_lic, jml_lic))      call nf95_inq_varid(ncid, 'longitude', varid)
233      ALLOCATE(dlon_lic(iml_lic))      call nf95_gw_var(ncid, varid, dlon_lic)
234      ALLOCATE(dlat_lic(jml_lic))      iml_lic = size(dlon_lic)
235    
236        call nf95_inq_varid(ncid, 'latitude', varid)
237        call nf95_gw_var(ncid, varid, dlat_lic)
238        jml_lic = size(dlat_lic)
239    
240        call nf95_inq_varid(ncid, 'landice', varid)
241      ALLOCATE(fraclic(iml_lic, jml_lic))      ALLOCATE(fraclic(iml_lic, jml_lic))
242      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)  
243    
244      ! Interpolation sur la grille T du modèle :      call nf95_close(ncid)
245      PRINT *, 'Dimensions de "landice"'  
246        ! Interpolation sur la grille T du mod\`ele :
247        PRINT *, 'Dimensions de "landiceref.nc"'
248      print *, "iml_lic = ", iml_lic      print *, "iml_lic = ", iml_lic
249      print *, "jml_lic = ", jml_lic      print *, "jml_lic = ", jml_lic
250    
251      ! Si les coordonnées sont en degrés, on les transforme :      ! Si les coordonn\'ees sont en degr\'es, on les transforme :
252      IF (MAXVAL( lon_lic(:, :) ) > pi)  THEN      IF (MAXVAL(dlon_lic) > pi) THEN
253         lon_lic(:, :) = lon_lic(:, :) * pi / 180.         dlon_lic = dlon_lic * pi / 180.
254      ENDIF      ENDIF
255      IF (maxval( lat_lic(:, :) ) > pi) THEN      IF (maxval(dlat_lic) > pi) THEN
256         lat_lic(:, :) = lat_lic(:, :) * pi/ 180.         dlat_lic = dlat_lic * pi/ 180.
257      ENDIF      ENDIF
258    
     dlon_lic(:) = lon_lic(:, 1)  
     dlat_lic(:) = lat_lic(1, :)  
   
259      flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &      flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
260           rlatu)           rlatu)
261      flic_tmp(iim + 1, :) = flic_tmp(1, :)      flic_tmp(iim + 1, :) = flic_tmp(1, :)
262    
263        deallocate(dlon_lic, dlat_lic) ! pointers
264    
265      ! Passage sur la grille physique      ! Passage sur la grille physique
266      pctsrf(:, :)=0.      pctsrf = 0.
267      pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)      pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
268      ! Adéquation avec le maque terre/mer      ! Ad\'equation avec le maque terre/mer
269      WHERE (pctsrf(:, is_lic) < EPSFRA ) pctsrf(:, is_lic) = 0.      WHERE (pctsrf(:, is_lic) < EPSFRA) pctsrf(:, is_lic) = 0.
270      WHERE (zmasq(:) < EPSFRA) pctsrf(:, is_lic) = 0.      WHERE (zmasq < EPSFRA) pctsrf(:, is_lic) = 0.
271      pctsrf(:, is_ter) = zmasq(:)      pctsrf(:, is_ter) = zmasq
272      where (zmasq(:) > EPSFRA)      where (zmasq > EPSFRA)
273         where (pctsrf(:, is_lic) >= zmasq(:))         where (pctsrf(:, is_lic) >= zmasq)
274            pctsrf(:, is_lic) = zmasq(:)            pctsrf(:, is_lic) = zmasq
275            pctsrf(:, is_ter) = 0.            pctsrf(:, is_ter) = 0.
276         elsewhere         elsewhere
277            pctsrf(:, is_ter) = zmasq(:) - pctsrf(:, is_lic)            pctsrf(:, is_ter) = zmasq - pctsrf(:, is_lic)
278            where (pctsrf(:, is_ter) < EPSFRA)            where (pctsrf(:, is_ter) < EPSFRA)
279               pctsrf(:, is_ter) = 0.               pctsrf(:, is_ter) = 0.
280               pctsrf(:, is_lic) = zmasq(:)               pctsrf(:, is_lic) = zmasq
281            end where            end where
282         end where         end where
283      end where      end where
284    
285      ! 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
286      ! de mer à 0) :      ! de mer \`a 0) :
287      pctsrf(:, is_oce) = 1. - zmasq(:)      pctsrf(:, is_oce) = 1. - zmasq
288      WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.      WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
289    
290      ! Vérification que somme des sous-surfaces vaut 1:      ! V\'erification que somme des sous-surfaces vaut 1 :
291      ji = count(abs(sum(pctsrf(:, :), dim = 2) - 1. ) > EPSFRA)      ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
292      IF (ji /= 0) PRINT *, 'Problème répartition sous maille pour', ji, 'points'      IF (ji /= 0) then
293           PRINT *, 'Bad surface percentages for ', ji, 'points'
294        end IF
295    
296      ! Calcul intermédiaire:      ! Calcul interm\'ediaire :
297      CALL massdair(p3d, masse)      CALL massdair(p3d, masse)
298    
299      print *, 'ALPHAX = ', alphax      print *, 'ALPHAX = ', alphax
300    
301      forall  (l = 1:llm)      forall (l = 1:llm)
302         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
303         masse(:, jjm + 1, l) = &         masse(:, jjm + 1, l) = &
304              SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols              SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
# Line 289  contains Line 306  contains
306    
307      ! Initialisation pour traceurs:      ! Initialisation pour traceurs:
308      call iniadvtrac      call iniadvtrac
     ! Ecriture:  
     CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &  
          tetagrot, tetatemp)  
     itau_dyn = 0  
309      itau_phy = 0      itau_phy = 0
310      day_ref = dayref      day_ref = dayref
311      annee_ref = anneeref      annee_ref = anneeref
312    
313      CALL geopot(ip1jmp1, tpot, pk , pks,  phis  , phi)      CALL geopot(teta, pk , pks, phis, phi)
314      CALL caldyn0(0, uvent, vvent, tpot, psol, masse, pk, phis, phi, w, &      CALL caldyn0(ucov, vcov, teta, ps, masse, pk, phis, phi, w, pbaru, &
315           pbaru, pbarv, 0)           pbarv)
316      CALL dynredem0("start.nc", dayref, phis)      CALL dynredem0("start.nc", dayref, phis)
317      CALL dynredem1("start.nc", 0., vvent, uvent, tpot, q3d, 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  
318    
319      ! Initialisations :      ! Initialisations :
320      tsolsrf(:, is_ter) = tsol      tsolsrf(:, is_ter) = tsol
# Line 326  contains Line 330  contains
330      albe(:, is_oce) = 0.5      albe(:, is_oce) = 0.5
331      albe(:, is_sic) = 0.6      albe(:, is_sic) = 0.6
332      alblw = albe      alblw = albe
333      evap(:, :) = 0.      evap = 0.
334      qsolsrf(:, is_ter) = 150.      qsolsrf(:, is_ter) = 150.
335      qsolsrf(:, is_lic) = 150.      qsolsrf(:, is_lic) = 150.
336      qsolsrf(:, is_oce) = 150.      qsolsrf(:, is_oce) = 150.
# Line 340  contains Line 344  contains
344      q_ancien = 0.      q_ancien = 0.
345      agesno = 0.      agesno = 0.
346      !IM "slab" ocean      !IM "slab" ocean
347      tslab(:) = tsolsrf(:, is_oce)      tslab = tsolsrf(:, is_oce)
348      seaice = 0.      seaice = 0.
349    
350      frugs(:, is_oce) = rugmer(:)      frugs(:, is_oce) = rugmer
351      frugs(:, is_ter) = MAX(1.e-05, zstd(:) * zsig(:) / 2)      frugs(:, is_ter) = MAX(1e-5, zstd * zsig / 2)
352      frugs(:, is_lic) = MAX(1.e-05, zstd(:) * zsig(:) / 2)      frugs(:, is_lic) = MAX(1e-5, zstd * zsig / 2)
353      frugs(:, is_sic) = 0.001      frugs(:, is_sic) = 0.001
354      fder = 0.      fder = 0.
355      clwcon = 0.      clwcon = 0.
356      rnebcon = 0.      rnebcon = 0.
357      ratqs = 0.      ratqs = 0.
358      run_off_lic_0 = 0.      run_off_lic_0 = 0.
359        sig1 = 0.
360        w01 = 0.
361    
362      call phyredem("startphy.nc", phystep, radpas, latfi, lonfi, pctsrf, &      call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
363           tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &           tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
364           evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &           evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
365           agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel, &           agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
366           t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)           t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
367      CALL histclo      CALL histclo
368    
369    END SUBROUTINE etat0    END SUBROUTINE etat0
Legend:
Removed from v.10  
changed lines
  Added in v.97
  ViewVC Help
Powered by ViewVC 1.1.21