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

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

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