source: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/etat0_netcdf.F @ 1160

Last change on this file since 1160 was 1154, checked in by lguez, 15 years ago

-- Added "NetCDF95" interface in "bibio".
-- NetCDF95 uses the module "typesizes", which is part of NetCDF, so we
exclude dependency on "typesizes" in "bld.cfg".
-- Added "assert_eq" and "assert" procedures, which are in the public
part of Numerical Recipes.
-- Added some interpolation and regridding utilities in "bibio".
-- Added the ability to read an ozone climatology from a NetCDF file.
-- Commented out unused variables and code in "etat0_netcdf".
-- Updated calls to NetCDF in "etat0_netcdf": from Fortran 77
interface to Fortran 90 interface.
-- Removed useless "deallocate" at the end of "etat0_netcdf".
-- Corrected some declarations not conforming to Fortran standard, such
as "integer*4", or obsolescent such as "character*4".
-- Replaced some calls to not-standard function "float" by calls to
"real".
-- On Brodie at IDRIS, the NetCDF library compiled with OpenMP should
be used. Changed path in "arch-SX8_BRODIE.path".
-- Added warning for incompatibility of debugging options and OpenMP
parallelization in "makelmdz_fcm".

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