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

Last change on this file since 155 was 115, checked in by lmdzadmin, 24 years ago

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