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

Diff of /trunk/dyn3d/etat0.f

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

trunk/libf/dyn3d/etat0.f90 revision 5 by guez, Mon Mar 3 16:32:04 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
47      use exner_hyb_m, only: exner_hyb      use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
48      use regr_coefoz_m, only: regr_coefoz      use regr_pr_o3_m, only: regr_pr_o3
49      use advtrac_m, only: iniadvtrac      use serre, only: alphax
50      use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf90_nowrite, &      use startdyn, only: start_init_dyn
51           nf90_get_var, handle_err      USE start_init_orog_m, only: start_init_orog, mask
52      use pressure_m, only: pls, p3d      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      integer ncid, varid, ncerr, month      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 131  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
206           call regr_pr_o3(p3d, q(:, :, :, 5))
207           ! Convert from mole fraction to mass fraction:
208           q(:, :, :, 5) = q(:, :, :, 5) * 48. / 29.
209        end if
210    
211         ! Compute ozone parameters on the LMDZ grid:      tsol = pack(tsol_2d, dyn_phy)
212         call regr_coefoz      qsol = pack(qsol_2d, dyn_phy)
213        sn = 0. ! snow
214        radsol = 0.
215        tslab = 0. ! IM "slab" ocean
216        seaice = 0.
217        rugmer = 0.001
218        zmea = pack(zmea_2d, dyn_phy)
219        zstd = pack(zstd_2d, dyn_phy)
220        zsig = pack(zsig_2d, dyn_phy)
221        zgam = pack(zgam_2d, dyn_phy)
222        zthe = pack(zthe_2d, dyn_phy)
223        zpic = pack(zpic_2d, dyn_phy)
224        zval = pack(zval_2d, dyn_phy)
225    
226         ! Find the month containing the day number "dayref":      ! On initialise les sous-surfaces.
227         month = (dayref - 1) / 30 + 1      ! Lecture du fichier glace de terre pour fixer la fraction de terre
228         print *, "month = ", month      ! et de glace de terre :
   
        call nf95_open("coefoz_LMDZ.nc", nf90_nowrite, ncid)  
   
        ! Get data at the right month from the input file:  
        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)  
   
        call nf95_close(ncid)  
        ! Latitudes are in increasing order in the input file while  
        ! "rlatu" is in decreasing order so we need to invert order. Also, we  
        ! compute mass fraction from mole fraction:  
        q3d(:, :, :, 5) = q3d(:, jjm+1:1:-1, :, 5)  * 48. / 29.  
     end if  
229    
230      tsol(:) = pack(tsol_2d, dyn_phy)      call nf95_open("landiceref.nc", nf90_nowrite, ncid)
     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)  
231    
232      rugsrel(:) = 0.      call nf95_inq_varid(ncid, 'longitude', varid)
233      IF (ok_orodr) rugsrel(:) = MAX(1.e-05, zstd(:) * zsig(:) / 2)      call nf95_gw_var(ncid, varid, dlon_lic)
234        iml_lic = size(dlon_lic)
235    
236      ! On initialise les sous-surfaces:      call nf95_inq_varid(ncid, 'latitude', varid)
237      ! Lecture du fichier glace de terre pour fixer la fraction de terre      call nf95_gw_var(ncid, varid, dlat_lic)
238      ! et de glace de terre:      jml_lic = size(dlat_lic)
239      CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, &  
240           ttm_tmp, fid)      call nf95_inq_varid(ncid, 'landice', varid)
     ALLOCATE(lat_lic(iml_lic, jml_lic))  
     ALLOCATE(lon_lic(iml_lic, jml_lic))  
     ALLOCATE(dlon_lic(iml_lic))  
     ALLOCATE(dlat_lic(jml_lic))  
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)
243           llm_tmp, lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt,  &  
244           fid)      call nf95_close(ncid)
     CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp &  
          , 1, 1, fraclic)  
     CALL flinclo(fid)  
245    
246      ! Interpolation sur la grille T du modèle :      ! Interpolation sur la grille T du mod\`ele :
247      PRINT *, 'Dimensions de "landice"'      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 307  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 344  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 358  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.5  
changed lines
  Added in v.97
  ViewVC Help
Powered by ViewVC 1.1.21