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

Last change on this file since 524 was 524, checked in by lmdzadmin, 20 years ago

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