source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/etat0_netcdf.F @ 174

Last change on this file since 174 was 174, checked in by lmdzadmin, 23 years ago

Modifs pour l'utilisation du masque venant de l'ocean en cas de couple
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.5 KB
Line 
1      SUBROUTINE etat0_netcdf
2   
3      USE startvar
4      USE ioipsl
5      !
6      IMPLICIT NONE
7      !
8#include "netcdf.inc"
9#include "dimensions.h"
10#include "paramet.h"
11      !
12      !
13c      INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2,
14c     .KLON=KFDIA-KIDIA+1,KLEV=llm
15      !
16#include "comgeom2.h"
17#include "comvert.h"
18#include "comconst.h"
19#include "indicesol.h"
20#include "dimphy.h"
21#include "dimsoil.h"
22      !
23      REAL :: latfi(klon), lonfi(klon)
24      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
25     . psol(iip1, jjp1), phis(iip1, jjp1)
26      REAL :: p3d(iip1, jjp1, llm+1)
27      REAL :: uvent(iip1, jjp1, llm)
28      REAL :: vvent(iip1, jjm, llm)
29      REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
30      REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
31      REAL :: tsol(klon), qsol(klon), sn(klon)
32      REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf)
33      REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
34      REAL :: tsoil(klon,nsoilmx,nbsrf)
35      REAL :: radsol(klon),rain_fall(klon), snow_fall(klon)
36      REAL :: solsw(klon), sollw(klon)
37      REAL :: deltat(klon), frugs(klon,nbsrf), agesno(klon),rugmer(klon)
38      REAL :: zmea(iip1*jjp1), zstd(iip1*jjp1)
39      REAL :: zsig(iip1*jjp1), zgam(iip1*jjp1), zthe(iip1*jjp1)
40      REAL :: zpic(iip1*jjp1), zval(iip1*jjp1), rugsrel(iip1*jjp1)
41      REAL :: qd(iip1, jjp1, llm)
42      REAL :: pctsrf(klon, nbsrf)
43      REAL :: t_ancien(klon,klev), q_ancien(klon,klev)      !
44      ! declarations pour lecture glace de mer
45      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
46      INTEGER :: itaul(1), fid
47      REAL :: lev(1), date, dt
48      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
49      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
50      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
51      REAL :: flic_tmp(iip1, jjp1)
52      REAL :: champint(iim, jjp1)
53      !
54
55      CHARACTER*80 :: varname
56      !
57      INTEGER :: i,j, ig, l, ji
58      REAL :: xpi
59      !
60      REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
61      REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
62      REAL :: workvar(iip1,jjp1,llm)
63      !
64      REAL ::  prefkap, unskap
65      !
66      REAL :: q_sat
67      EXTERNAL q_sat
68      real :: time_step,t_ops,t_wrt
69
70#include "comdissnew.h"
71#include "control.h"
72#include "serre.h"
73#include "clesph0.h"
74
75      INTEGER  ::        longcles
76      PARAMETER      ( longcles  = 20 )
77      REAL :: clesphy0 ( longcles       )
78      REAL :: p(iip1,jjp1,llm)
79      INTEGER :: itau, iday
80      REAL :: masse(iip1,jjp1,llm)
81      REAL :: xpn,xps,xppn(iim),xpps(iim)
82      real :: time
83      REAL :: phi(ip1jmp1,llm)
84      REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
85      REAL :: w(ip1jmp1,llm)
86      REAL ::phystep,co2_ppm,solaire
87      INTEGER :: radpas
88
89      CHARACTER*80 :: visu_file
90      INTEGER :: visuid
91
92! pour la lecture du fichier masque ocean
93      integer :: nid_o2a
94      logical :: couple = .false.
95      INTEGER :: iml_omask, jml_omask
96      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
97      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
98      REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
99      real, dimension(klon) :: ocemask_fi
100      integer :: isst(klon-2)
101
102      !
103      !   Constantes
104      !
105      pi     = 4. * ATAN(1.)
106      rad    = 6371229.
107      omeg   = 4.* ASIN(1.)/(24.*3600.)
108      g      = 9.8
109      daysec = 86400.
110      kappa  = 0.2857143
111      cpp    = 1004.70885
112      !
113      preff     = 101325.
114      unskap = 1./kappa
115      !
116      jmp1    = jjm + 1
117      !
118      !    Construct a grid
119      !
120
121      CALL defrun_new(99,.TRUE.,clesphy0)
122
123      dtvr   = daysec/FLOAT(day_step)
124      print*,'dtvr',dtvr
125
126      CALL inicons0()
127      CALL inigeom()
128      !
129      CALL inifilr()
130      !
131      latfi(1) = ASIN(1.0)
132      DO j = 2, jjm
133        DO i = 1, iim
134          latfi((j-2)*iim+1+i)=  rlatu(j)
135        ENDDO
136      ENDDO
137      latfi(klon) = - ASIN(1.0)
138      !
139      lonfi(1) = 0.0
140      DO j = 2, jjm
141        DO i = 1, iim
142          lonfi((j-2)*iim+1+i) =  rlonv(i)
143        ENDDO
144      ENDDO
145      lonfi(klon) = 0.0
146      !
147      xpi = 2.0 * ASIN(1.0)
148      DO ig = 1, klon
149        latfi(ig) = latfi(ig) * 180.0 / xpi
150        lonfi(ig) = lonfi(ig) * 180.0 / xpi
151      ENDDO
152      !
153      varname = 'relief'
154      ! This line needs to be replaced by a call to restget to get the values in the restart file
155      orog(:,:) = 0.0
156       CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0)
157      !
158      WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
159      WRITE(*,'(49I1)') INT(orog(:,:))
160      !
161      varname = 'rugosite'
162      ! This line needs to be replaced by a call to restget to get the values in the restart file
163      rugo(:,:) = 0.0
164       CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0)
165      !
166      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite'
167      WRITE(*,'(49I1)') INT(rugo(:,:)*10)
168      !
169      varname = 'masque'
170      ! This line needs to be replaced by a call to restget to get the values in the restart file
171      masque(:,:) = 0.0
172       CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0)
173      !
174      WRITE(*,*) 'MASQUE construit : Masque'
175      WRITE(*,'(49I1)') INT(masque(:,:))
176      !
177      !
178
179
180C
181C on initialise les sous surfaces
182C
183      pctsrf=0.
184      !cree le masque a partir du fichier relief
185      varname = 'zmasq'
186      zmasq(:) = 0.
187      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0)
188      WHERE (zmasq(1 : klon) .LE. EPSFRA)
189          zmasq(1 : klon) = 0.
190      END WHERE
191      WRITE(*,*)zmasq
192
193
194
195
196      varname = 'psol'
197      psol(:,:) = 0.0
198      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0)
199      !
200      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
201      !  anyway.
202      !
203      WRITE(*,*) 'PSOL :', psol(10,20)
204      WRITE(*,*) ap(:), bp(:)
205      CALL pression(ip1jmp1, ap, bp, psol, p3d)
206      WRITE(*,*) 'P3D :', p3d(10,20,:)
207      CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
208      WRITE(*,*) 'PK:', pk(10,20,:)
209      !
210      !
211      !
212      prefkap =  preff  ** kappa
213      WRITE(*,*) 'unskap, cpp,  preff :', unskap, cpp,  preff
214      DO l = 1, llm
215        DO j=1,jjp1
216          DO i =1, iip1
217            pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
218           ENDDO
219        ENDDO
220      ENDDO
221      !
222      WRITE(*,*) 'PLS :', pls(10,20,:)
223      !
224      varname = 'surfgeo'
225      phis(:,:) = 0.0
226      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0)
227      write(*,*) 'Phis = '
228      write(*,*)phis
229      !
230      varname = 'u'
231      uvent(:,:,:) = 0.0
232      CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
233     . workvar, uvent, 0.0)
234      ! 
235      varname = 'v'
236      vvent(:,:,:) = 0.0
237      CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
238     . workvar, vvent, 0.0)
239      !
240      varname = 't'
241      t3d(:,:,:) = 0.0
242      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
243     . workvar, t3d, 0.0)
244      !
245      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
246     .                          maxval(t3d(:,:,:))
247      varname = 'tpot'
248      tpot(:,:,:) = 0.0
249      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
250     . pk, tpot, 0.0)
251      !
252      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
253     .                          maxval(t3d(:,:,:))
254      WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)),
255     .                          maxval(pls(:,:,:))
256      DO l = 1, llm
257        DO j=1,jjp1
258          DO i =1, iip1-1
259           qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
260          ENDDO
261          qsat(iip1,j,l) = qsat(1,j,l)
262        ENDDO
263      ENDDO
264      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
265     .                           maxval(qsat(:,:,:))
266      !
267      WRITE(*,*) 'QSAT :', qsat(10,20,:)
268      !
269      varname = 'q'
270      q3d(:,:,:,:) = 0.0
271      qd(:,:,:) = 0.0
272      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
273     .                           maxval(qsat(:,:,:))
274      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
275     . qsat, qd, 0.0)
276      q3d(:,:,:,1) = qd(:,:,:)
277      !
278      varname = 'tsol'
279      ! This line needs to be replaced by a call to restget to get the values in the restart file
280      tsol(:) = 0.0
281      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol,0.0)
282      !
283      WRITE(*,*) 'TSOL construit :'
284      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
285      !
286      varname = 'qsol'
287      qsol(:) = 0.0
288      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol,0.0)
289      !
290      varname = 'snow'
291      sn(:) = 0.0
292      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn,0.0)
293      !
294      varname = 'rads'
295      radsol(:) = 0.0
296      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0)
297      !
298      varname = 'deltat'
299      deltat(:) = 0.0
300      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0)
301      !
302      varname = 'rugmer'
303      rugmer(:) = 0.0
304      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0)
305      !
306      varname = 'agsno'
307      agesno(:) = 0.0
308      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0)
309
310      varname = 'zmea'
311      zmea(:) = 0.0
312      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0)
313      varname = 'zstd'
314      zstd(:) = 0.0
315      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0)
316      varname = 'zsig'
317      zsig(:) = 0.0
318      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0)
319      varname = 'zgam'
320      zgam(:) = 0.0
321      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0)
322      varname = 'zthe'
323      zthe(:) = 0.0
324      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0)
325      varname = 'zpic'
326      zpic(:) = 0.0
327      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0)
328      varname = 'zval'
329      zval(:) = 0.0
330      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0)
331      rugsrel(:) = 0.0
332
333C
334C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
335C utilise pour calculer les poids et pour assurer l'adequation entre les
336C fractions d'ocean vu par l'atmosphere et l'ocean
337C
338
339      write(*,*)'Essai de lecture masque ocean'
340      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
341      if (iret .ne. 0) then
342        write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
343        write(*,*)'Run force'
344      else
345        couple = .true.
346        iret = nf_close(nid_o2a)
347        call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
348     $    , nid_o2a)
349        if (iml_omask /= iim .or. jml_omask /= jjp1) then
350          write(*,*)'Dimensions non compatibles pour masque ocean'
351          write(*,*)'iim = ',iim,' iml_omask = ',iml_omask
352          write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
353          stop
354        endif
355        ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
356        ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
357        ALLOCATE(dlon_omask(iml_omask), stat=iret)
358        ALLOCATE(dlat_omask(jml_omask), stat=iret)
359        ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
360        ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
361        CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
362     $    , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
363        CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
364     $      ttm_tmp, 1, 1, ocetmp)
365        CALL flinclo(fid)
366        dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
367        dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
368        ocemask = ocetmp
369        if (dlat_omask(1) < dlat_omask(jml_omask)) then
370          do j = 1, jml_omask
371            ocemask(:,j) = ocetmp(:,jml_omask-j+1)
372          enddo
373        endif
374C
375C passage masque ocean a la grille physique
376C
377        ocemask_fi(1) = ocemask(1,1)
378        do j = 2, jjm
379          do i = 1, iim
380            ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
381          enddo
382        enddo
383        ocemask_fi(klon) = ocemask(1,jjp1)
384        zmasq = 1. - ocemask_fi
385        isst = 0
386        where (ocemask_fi(2:klon-1) >0.) isst = 1
387        write(45,'(72i1)')isst
388      endif
389
390
391C
392C lecture du fichier glace de terre pour fixer la fraction de terre
393C et de glace de terre
394C
395      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
396     $    , fid)
397      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
398      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
399      ALLOCATE(dlon_lic(iml_lic), stat=iret)
400      ALLOCATE(dlat_lic(jml_lic), stat=iret)
401      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
402      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
403     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
404      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
405     $    , 1, 1, fraclic)
406      CALL flinclo(fid)
407C
408C interpolation sur la grille T du modele
409C
410      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ',
411     $    iml_lic, jml_lic
412c
413C sil les coordonnees sont en degres, on les transforme
414C
415      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
416          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
417      ENDIF
418      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
419          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
420      ENDIF
421
422      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
423      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
424C
425      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
426     $    ,iim, jjp1,
427     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
428c$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
429      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
430C
431C passage sur la grille physique
432C
433      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
434     $    pctsrf(1:klon, is_lic))
435C adequation avec le maque terre/mer
436      WHERE (pctsrf(1 : klon, is_lic) .LE. EPSFRA )
437          pctsrf(1 : klon, is_lic) = 0.
438      END WHERE
439      WHERE (zmasq( 1 : klon) .LE. EPSFRA)
440          pctsrf(1 : klon, is_lic) = 0.
441      END WHERE
442      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
443      DO ji = 1, klon
444        IF (zmasq(ji) .GT. EPSFRA) THEN
445            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
446                pctsrf(ji, is_lic) = zmasq(ji)
447                pctsrf(ji, is_ter) = 0.
448            ELSE
449                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
450            ENDIF
451        ENDIF
452      END DO
453C
454C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
455C
456      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
457
458
459      WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
460          pctsrf(1 : klon, is_oce) = 0.
461      END WHERE
462
463      if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
464
465      isst = 0
466      where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
467      write(46,'(72i1)')isst
468C
469C verif que somme des sous surface = 1
470C
471      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
472     $    .GT. EPSFRA)
473      IF (ji .NE. 0) THEN
474          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
475      ENDIF
476
477
478
479
480
481C Calcul intermediaire
482c
483      CALL massdair( p3d, masse  )
484c
485
486      print *,' ALPHAX ',alphax
487
488      DO  l = 1, llm
489        DO  i    = 1, iim
490          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
491          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
492        ENDDO
493          xpn      = SUM(xppn)/apoln
494          xps      = SUM(xpps)/apols
495        DO i   = 1, iip1
496          masse(   i   ,   1     ,  l )   = xpn
497          masse(   i   ,   jjp1  ,  l )   = xps
498        ENDDO
499      ENDDO
500      q3d(iip1,:,:,:) = q3d(1,:,:,:)
501      phis(iip1,:) = phis(1,:)
502
503C Ecriture
504
505
506      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
507     *                tetagdiv, tetagrot , tetatemp              )
508      print*,'sortie inidissip'
509      itau = 0
510      iday = dayref +itau/day_step
511      time = FLOAT(itau-(iday-dayref)*day_step)/day_step
512c     
513      IF(time.GT.1)  THEN
514       time = time - 1
515       iday = iday + 1
516      ENDIF
517      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
518      print*,'sortie geopot'
519     
520      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
521     *                phi,w, pbaru,pbarv,time+iday-dayref   )
522       print*,'sortie caldyn0'     
523      CALL dynredem0("start.nc",dayref,anneeref,phis,nqmx)
524      print*,'sortie dynredem0'
525      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
526     .                            psol)
527      print*,'sortie dynredem1'
528C
529C Ecriture etat initial physique
530C
531      phystep   = dtvr * FLOAT(iphysiq)
532      radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
533      co2_ppm   = 330.0
534      solaire   = 1370.0
535
536c      call physdem(lonfi, latfi, phystep,radpas,co2_ppm,
537c     .                   solaire,tsol, qsol,
538c     .                   sn, radsol, deltat, rugmer,
539c     .                   agesno, zmea, zstd, zsig,
540c     .                   zgam, zthe, zpic, zval,
541c     .                   rugsrel)
542
543c
544c Initialisation
545c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
546c
547      tsolsrf(:,is_ter) = tsol
548      tsolsrf(:,is_lic) = tsol
549      tsolsrf(:,is_oce) = tsol
550      tsolsrf(:,is_sic) = tsol
551      snsrf(:,is_ter) = sn
552      snsrf(:,is_lic) = sn
553      snsrf(:,is_oce) = sn
554      snsrf(:,is_sic) = sn
555      albe(:,is_ter) = 0.08
556      albe(:,is_lic) = 0.6
557      albe(:,is_oce) = 0.5
558      albe(:,is_sic) = 0.6
559      evap(:,:) = 0.
560      qsolsrf(:,is_ter) = qsol
561      qsolsrf(:,is_lic) = qsol
562      qsolsrf(:,is_oce) = 150.
563      qsolsrf(:,is_sic) = 150.
564      do i = 1, nbsrf
565        do j = 1, nsoilmx
566          tsoil(:,j,i) = tsol
567        enddo
568      enddo
569      rain_fall = 0.; snow_fall = 0.
570      solsw = 165.
571      sollw = -53.
572      t_ancien = 273.15
573      q_ancien = 0.
574      agesno = 0.
575      deltat = 0.
576      frugs(:,is_oce) = rugmer
577      frugs(:,is_ter) = rugmer
578      frugs(:,is_lic) = rugmer
579      frugs(:,is_sic) = rugmer
580
581      call physdem("startphy.nc",phystep,radpas, co2_ppm, solaire,
582     $    latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, snsrf,
583     $    albe, evap, rain_fall, snow_fall, solsw, sollw,
584     $    radsol, frugs,  agesno,
585     $    zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
586     $    t_ancien, q_ancien)
587
588C     Sortie Visu pour les champs dynamiques
589      print*,'sortie visu'
590      time_step = 1.
591      t_ops = 2.
592      t_wrt = 2.
593      itau = 2.
594      visu_file='Etat0_visu.nc'
595      CALL initdynav(visu_file,dayref,anneeref,time_step,
596     .              t_ops, t_wrt, nqmx, visuid)
597      CALL writedynav(visuid, nqmx, itau,vvent ,
598     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
599      print*,'entree histclo'
600      CALL histclo
601      RETURN
602      !
603      END SUBROUTINE etat0_netcdf
Note: See TracBrowser for help on using the repository browser.