source: LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F @ 762

Last change on this file since 762 was 762, checked in by Laurent Fairhead, 17 years ago

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