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

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

Lots of stuff, plus particulierement:

  • appel a ORCHIDEE en etat de marche (pb de grille subsiste)
  • modifs de Pascale sur soil dans le cas ou ok_veget=false -
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.7 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) .LT. EPSFRA)
189          zmasq(1 : klon) = 0.
190      END WHERE
191      WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
192          zmasq(1 : klon) = 1.
193      END WHERE
194      WRITE(*,*)zmasq
195
196
197
198
199      varname = 'psol'
200      psol(:,:) = 0.0
201      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0)
202      !
203      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
204      !  anyway.
205      !
206      WRITE(*,*) 'PSOL :', psol(10,20)
207      WRITE(*,*) ap(:), bp(:)
208      CALL pression(ip1jmp1, ap, bp, psol, p3d)
209      WRITE(*,*) 'P3D :', p3d(10,20,:)
210      CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
211      WRITE(*,*) 'PK:', pk(10,20,:)
212      !
213      !
214      !
215      prefkap =  preff  ** kappa
216      WRITE(*,*) 'unskap, cpp,  preff :', unskap, cpp,  preff
217      DO l = 1, llm
218        DO j=1,jjp1
219          DO i =1, iip1
220            pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
221           ENDDO
222        ENDDO
223      ENDDO
224      !
225      WRITE(*,*) 'PLS :', pls(10,20,:)
226      !
227      varname = 'surfgeo'
228      phis(:,:) = 0.0
229      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0)
230      write(*,*) 'Phis = '
231      write(*,*)phis
232      !
233      varname = 'u'
234      uvent(:,:,:) = 0.0
235      CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
236     . workvar, uvent, 0.0)
237      ! 
238      varname = 'v'
239      vvent(:,:,:) = 0.0
240      CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
241     . workvar, vvent, 0.0)
242      !
243      varname = 't'
244      t3d(:,:,:) = 0.0
245      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
246     . workvar, t3d, 0.0)
247      !
248      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
249     .                          maxval(t3d(:,:,:))
250      varname = 'tpot'
251      tpot(:,:,:) = 0.0
252      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
253     . pk, tpot, 0.0)
254      !
255      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
256     .                          maxval(t3d(:,:,:))
257      WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)),
258     .                          maxval(pls(:,:,:))
259      DO l = 1, llm
260        DO j=1,jjp1
261          DO i =1, iip1-1
262           qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
263          ENDDO
264          qsat(iip1,j,l) = qsat(1,j,l)
265        ENDDO
266      ENDDO
267      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
268     .                           maxval(qsat(:,:,:))
269      !
270      WRITE(*,*) 'QSAT :', qsat(10,20,:)
271      !
272      varname = 'q'
273      q3d(:,:,:,:) = 0.0
274      qd(:,:,:) = 0.0
275      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
276     .                           maxval(qsat(:,:,:))
277      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
278     . qsat, qd, 0.0)
279      q3d(:,:,:,1) = qd(:,:,:)
280      !
281      varname = 'tsol'
282      ! This line needs to be replaced by a call to restget to get the values in the restart file
283      tsol(:) = 0.0
284      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol,0.0)
285      !
286      WRITE(*,*) 'TSOL construit :'
287      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
288      !
289      varname = 'qsol'
290      qsol(:) = 0.0
291      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol,0.0)
292      !
293      varname = 'snow'
294      sn(:) = 0.0
295      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn,0.0)
296      !
297      varname = 'rads'
298      radsol(:) = 0.0
299      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0)
300      !
301      varname = 'deltat'
302      deltat(:) = 0.0
303      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0)
304      !
305      varname = 'rugmer'
306      rugmer(:) = 0.0
307      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0)
308      !
309      varname = 'agsno'
310      agesno(:) = 0.0
311      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0)
312
313      varname = 'zmea'
314      zmea(:) = 0.0
315      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0)
316      varname = 'zstd'
317      zstd(:) = 0.0
318      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0)
319      varname = 'zsig'
320      zsig(:) = 0.0
321      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0)
322      varname = 'zgam'
323      zgam(:) = 0.0
324      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0)
325      varname = 'zthe'
326      zthe(:) = 0.0
327      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0)
328      varname = 'zpic'
329      zpic(:) = 0.0
330      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0)
331      varname = 'zval'
332      zval(:) = 0.0
333      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0)
334      rugsrel(:) = 0.0
335
336C
337C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
338C utilise pour calculer les poids et pour assurer l'adequation entre les
339C fractions d'ocean vu par l'atmosphere et l'ocean
340C
341
342      write(*,*)'Essai de lecture masque ocean'
343      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
344      if (iret .ne. 0) then
345        write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
346        write(*,*)'Run force'
347      else
348        couple = .true.
349        iret = nf_close(nid_o2a)
350        call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
351     $    , nid_o2a)
352        if (iml_omask /= iim .or. jml_omask /= jjp1) then
353          write(*,*)'Dimensions non compatibles pour masque ocean'
354          write(*,*)'iim = ',iim,' iml_omask = ',iml_omask
355          write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
356          stop
357        endif
358        ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
359        ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
360        ALLOCATE(dlon_omask(iml_omask), stat=iret)
361        ALLOCATE(dlat_omask(jml_omask), stat=iret)
362        ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
363        ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
364        CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
365     $    , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
366        CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
367     $      ttm_tmp, 1, 1, ocetmp)
368        CALL flinclo(fid)
369        dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
370        dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
371        ocemask = ocetmp
372        if (dlat_omask(1) < dlat_omask(jml_omask)) then
373          do j = 1, jml_omask
374            ocemask(:,j) = ocetmp(:,jml_omask-j+1)
375          enddo
376        endif
377C
378C passage masque ocean a la grille physique
379C
380        ocemask_fi(1) = ocemask(1,1)
381        do j = 2, jjm
382          do i = 1, iim
383            ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
384          enddo
385        enddo
386        ocemask_fi(klon) = ocemask(1,jjp1)
387        zmasq = 1. - ocemask_fi
388        isst = 0
389        where (ocemask_fi(2:klon-1) >0.) isst = 1
390        write(45,'(72i1)')isst
391      endif
392
393
394C
395C lecture du fichier glace de terre pour fixer la fraction de terre
396C et de glace de terre
397C
398      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
399     $    , fid)
400      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
401      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
402      ALLOCATE(dlon_lic(iml_lic), stat=iret)
403      ALLOCATE(dlat_lic(jml_lic), stat=iret)
404      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
405      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
406     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
407      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
408     $    , 1, 1, fraclic)
409      CALL flinclo(fid)
410C
411C interpolation sur la grille T du modele
412C
413      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ',
414     $    iml_lic, jml_lic
415c
416C sil les coordonnees sont en degres, on les transforme
417C
418      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
419          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
420      ENDIF
421      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
422          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
423      ENDIF
424
425      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
426      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
427C
428      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
429     $    ,iim, jjp1,
430     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
431c$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
432      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
433C
434C passage sur la grille physique
435C
436      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
437     $    pctsrf(1:klon, is_lic))
438C adequation avec le maque terre/mer
439c      zmasq(157) = 0.
440      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
441          pctsrf(1 : klon, is_lic) = 0.
442      END WHERE
443      WHERE (zmasq( 1 : klon) .LT. EPSFRA)
444          pctsrf(1 : klon, is_lic) = 0.
445      END WHERE
446      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
447      DO ji = 1, klon
448        IF (zmasq(ji) .GT. EPSFRA) THEN
449            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
450                pctsrf(ji, is_lic) = zmasq(ji)
451                pctsrf(ji, is_ter) = 0.
452            ELSE
453                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
454                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
455                    pctsrf(ji,is_ter) = 0.
456                    pctsrf(ji, is_lic) = zmasq(ji)
457                ENDIF
458            ENDIF
459        ENDIF
460      END DO
461C
462C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
463C
464      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
465
466
467      WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
468          pctsrf(1 : klon, is_oce) = 0.
469      END WHERE
470
471      if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
472
473      isst = 0
474      where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
475      write(46,'(72i1)')isst
476C
477C verif que somme des sous surface = 1
478C
479      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
480     $    .GT. EPSFRA)
481      IF (ji .NE. 0) THEN
482          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
483      ENDIF
484
485
486
487
488
489C Calcul intermediaire
490c
491      CALL massdair( p3d, masse  )
492c
493
494      print *,' ALPHAX ',alphax
495
496      DO  l = 1, llm
497        DO  i    = 1, iim
498          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
499          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
500        ENDDO
501          xpn      = SUM(xppn)/apoln
502          xps      = SUM(xpps)/apols
503        DO i   = 1, iip1
504          masse(   i   ,   1     ,  l )   = xpn
505          masse(   i   ,   jjp1  ,  l )   = xps
506        ENDDO
507      ENDDO
508      q3d(iip1,:,:,:) = q3d(1,:,:,:)
509      phis(iip1,:) = phis(1,:)
510
511C Ecriture
512
513
514      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
515     *                tetagdiv, tetagrot , tetatemp              )
516      print*,'sortie inidissip'
517      itau = 0
518      iday = dayref +itau/day_step
519      time = FLOAT(itau-(iday-dayref)*day_step)/day_step
520c     
521      IF(time.GT.1)  THEN
522       time = time - 1
523       iday = iday + 1
524      ENDIF
525      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
526      print*,'sortie geopot'
527     
528      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
529     *                phi,w, pbaru,pbarv,time+iday-dayref   )
530       print*,'sortie caldyn0'     
531      CALL dynredem0("start.nc",dayref,anneeref,phis,nqmx)
532      print*,'sortie dynredem0'
533      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
534     .                            psol)
535      print*,'sortie dynredem1'
536C
537C Ecriture etat initial physique
538C
539      phystep   = dtvr * FLOAT(iphysiq)
540      radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
541      co2_ppm   = 330.0
542      solaire   = 1370.0
543
544c      call physdem(lonfi, latfi, phystep,radpas,co2_ppm,
545c     .                   solaire,tsol, qsol,
546c     .                   sn, radsol, deltat, rugmer,
547c     .                   agesno, zmea, zstd, zsig,
548c     .                   zgam, zthe, zpic, zval,
549c     .                   rugsrel)
550
551c
552c Initialisation
553c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
554c
555      tsolsrf(:,is_ter) = tsol
556      tsolsrf(:,is_lic) = tsol
557      tsolsrf(:,is_oce) = tsol
558      tsolsrf(:,is_sic) = tsol
559      snsrf(:,is_ter) = sn
560      snsrf(:,is_lic) = sn
561      snsrf(:,is_oce) = sn
562      snsrf(:,is_sic) = sn
563      albe(:,is_ter) = 0.08
564      albe(:,is_lic) = 0.6
565      albe(:,is_oce) = 0.5
566      albe(:,is_sic) = 0.6
567      evap(:,:) = 0.
568      qsolsrf(:,is_ter) = 150
569      qsolsrf(:,is_lic) = 150
570      qsolsrf(:,is_oce) = 150.
571      qsolsrf(:,is_sic) = 150.
572      do i = 1, nbsrf
573        do j = 1, nsoilmx
574          tsoil(:,j,i) = tsol
575        enddo
576      enddo
577      rain_fall = 0.; snow_fall = 0.
578      solsw = 165.
579      sollw = -53.
580      t_ancien = 273.15
581      q_ancien = 0.
582      agesno = 0.
583      deltat = 0.
584      frugs(:,is_oce) = rugmer
585      frugs(:,is_ter) = rugmer
586      frugs(:,is_lic) = rugmer
587      frugs(:,is_sic) = 0.001
588
589      call physdem("startphy.nc",phystep,radpas, co2_ppm, solaire,
590     $    latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, snsrf,
591     $    albe, evap, rain_fall, snow_fall, solsw, sollw,
592     $    radsol, frugs,  agesno,
593     $    zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
594     $    t_ancien, q_ancien)
595
596C     Sortie Visu pour les champs dynamiques
597      print*,'sortie visu'
598      time_step = 1.
599      t_ops = 2.
600      t_wrt = 2.
601      itau = 2.
602      visu_file='Etat0_visu.nc'
603      CALL initdynav(visu_file,dayref,anneeref,time_step,
604     .              t_ops, t_wrt, nqmx, visuid)
605      CALL writedynav(visuid, nqmx, itau,vvent ,
606     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
607      print*,'entree histclo'
608      CALL histclo
609      RETURN
610      !
611      END SUBROUTINE etat0_netcdf
Note: See TracBrowser for help on using the repository browser.