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

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

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