source: LMDZ.3.3/trunk/libf/phylmd/1DUTILS.h @ 346

Last change on this file since 346 was 248, checked in by lmdz, 23 years ago

Routines pour le 1D F.Hourdin
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.6 KB
Line 
1      SUBROUTINE initial0(n,x)
2      IMPLICIT NONE
3      INTEGER n,i
4      REAL x(n)
5      DO 10 i=1,n
6         x(i)=0.
710    CONTINUE
8      RETURN
9      END
10      SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
11      IMPLICIT NONE
12c=======================================================================
13c   passage d'un champ de la grille scalaire a la grille physique
14c=======================================================================
15
16c-----------------------------------------------------------------------
17c   declarations:
18c   -------------
19
20      INTEGER im,jm,ngrid,nfield
21      REAL pdyn(im,jm,nfield)
22      REAL pfi(ngrid,nfield)
23
24      INTEGER i,j,ifield,ig
25
26c-----------------------------------------------------------------------
27c   calcul:
28c   -------
29
30      DO ifield=1,nfield
31c   traitement des poles
32         DO i=1,im
33            pdyn(i,1,ifield)=pfi(1,ifield)
34            pdyn(i,jm,ifield)=pfi(ngrid,ifield)
35         ENDDO
36
37c   traitement des point normaux
38         DO j=2,jm-1
39            ig=2+(j-2)*(im-1)
40            CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
41            pdyn(im,j,ifield)=pdyn(1,j,ifield)
42         ENDDO
43      ENDDO
44
45      RETURN
46      END
47
48      SUBROUTINE disvert1d(llm,kappa,sig,dsig,s,ds,dsig1,sdsig)
49      IMPLICIT NONE
50c
51c=======================================================================
52c
53c
54c    s = sigma ** kappa   :  coordonnee  verticale
55c    dsig(l)            : epaisseur de la couche l ds la coord.  s
56c    sig(l)             : sigma a l'interface des couches l et l-1
57c    ds(l)              : distance entre les couches l et l-1 en coord.s
58c
59c=======================================================================
60c
61c   declarations:
62c   -------------
63c
64      integer llm
65      real kappa,pi,x
66      real sig(llm+1),dsig(llm),s(llm),ds(llm),dsig1(llm),sdsig(llm)
67c
68      integer ll,l,lllm,lllmm1,lllmp1
69      real abid,abid2,som,quoi,quand,snorm,sigbid,sbid
70      REAL alpha,beta,h,zd,dz0,dz1
71      REAL gama,delta,deltaz,np
72
73      real nhaut
74      INTEGER ierr,ierr1,ierr2
75      real puiss
76
77      real asig,bsig,csig,esig,zsig,p,zz,sig1
78      REAL z1,z2
79c
80c-----------------------------------------------------------------------
81c
82      lllm=llm
83      lllmm1=lllm-1
84      lllmp1=lllm+1
85      pi=2.*asin(1.)
86
87      OPEN(99,file='sigma.def',status='old',form='formatted',
88     s   iostat=ierr1)
89      if(ierr1.ne.0) then
90         close(99)
91         open(99,file='esasig.def',status='old',form='formatted',
92     s   iostat=ierr2)
93      endif
94
95
96c-----------------------------------------------------------------------
97c   cas 1 on lit les options dans sigma.def:
98c   ----------------------------------------
99
100
101      if (ierr1.eq.0) then
102         PRINT*,'WARNING Lecture de sigma.def'
103         READ(99,*) deltaz
104         READ(99,*) h
105         READ(99,*) beta
106         READ(99,*) gama
107         READ(99,*) delta
108         READ(99,*) np
109         CLOSE(99)
110         alpha=deltaz/(llm*h)
111          do l= 1, llm
112             dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*
113     $          ( (tanh(gama*l)/tanh(gama*llm))**np +
114     $            (1.-l/FLOAT(llm))*delta )
115          enddo
116
117          sig(1)=1.
118          do l=1,llm-1
119             sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
120          enddo
121          sig(llm+1)=0.
122
123          do l = 1, llm
124             dsig(l) = sig(l)-sig(l+1)
125          enddo
126
127      else if(ierr2.eq.0) then
128
129         PRINT*,'WARNING Lecture de esasig.def'
130         READ(99,*) h
131         READ(99,*) dz0
132         READ(99,*) dz1
133         READ(99,*) nhaut
134         CLOSE(99)
135
136         dz0=dz0/h
137         dz1=dz1/h
138
139         sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
140
141         esig=1.
142
143         PRINT*
144         do l=1,20
145            print*,'esig=',esig
146            esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.)
147         enddo
148         PRINT*
149         csig=(1./sig1-1.)/(exp(esig)-1.)
150
151         DO L = 2, llm
152           zz=csig*(exp(esig*(l-1.))-1.)
153           sig(l) =1./(1.+zz)
154     &    * tanh(.5*(llm+1-l)/nhaut)
155         ENDDO
156         sig(1)=1.
157         sig(llm+1)=0.
158
159         do  l = 1, llm
160         dsig(l) =sig(l)-sig(l+1)
161         enddo
162
163      else
164
165         print*,'WARNING!!! Ancienne discretisation verticale'
166c        stop
167         h=7.
168         snorm  = 0.
169         do l = 1, llm
170            x = 2.*asin(1.) * (float(l)-0.5) / float(llm+1)
171            dsig(l) = 1.0 + 7.0 * sin(x)**2
172            snorm = snorm + dsig(l)
173         enddo
174         snorm = 1./snorm
175         do l = 1, llm
176            dsig(l) = dsig(l)*snorm
177         enddo
178         sig(llm+1) = 0.
179         do l = llm, 1, -1
180            sig(l) = sig(l+1) + dsig(l)
181         enddo
182
183      endif
184
185c-----------------------------------------------------------------------
186c   calcul de s, ds, sdsig...
187c   -------------------------
188
189       quoi      = 1. + 2.* kappa
190       s( llm )  = 1.
191       s(lllmm1) = quoi
192       IF( llm.gt.2 )  THEN
193          DO  ll = 2, lllmm1
194             l         = lllmp1 - ll
195             quand     = sig(l+1)/ sig(l)
196             s(l-1)    = quoi * (1.-quand) * s(l)  + quand * s(l+1)
197          ENDDO
198       END IF
199c
200       snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2)
201       DO l = 1, llm
202          s(l)    = s(l)/ snorm
203       ENDDO
204
205       DO l = 2, llm
206          ds(l)   = s(l-1) - s(l)
207       ENDDO
208       ds(1)  = 1. - s(1)
209c
210       DO l  = 1, llm
211          sdsig(l) = s(l) * dsig(l)
212          dsig1(l)= 1./dsig(l)
213       ENDDO
214
215c-----------------------------------------------------------------------
216c
217c     Diagnostique sur la discretisation verticale:
218c     ---------------------------------------------
219c
220      print*,'Diagnostique de la discretisation verticale'
221      print*
222      print*,'comparaison de sig(l) et (s(l)+s(l+1))/2)**(1/K)'
223      do 14 l=1,llm-1
224         sigbid=(0.5*(s(l)+s(l+1)))**(1./kappa)
225         print*,'sig(',l+1,')  = ',sig(l+1),
226     S           '    valeur approchee :',sigbid,'   ',dsig(l)
22714    continue
228      print*
229      print*,'comparaison de s(l) et (sig(l)+sig(l+1))/2)**K'
230      do 15 l=1,llm
231         sbid=(0.5*(sig(l+1)+sig(l)))**kappa
232         print*,'  s(',l,')  = ',s(l),
233     S           '    valeur approchee :',sbid
23415    continue
235c
236      PRINT*,'Altitude approchee z,dz'
237      PRINT*
238      z1=0.
239      print*,'   l       Z      DZ      Ztop   dsig'
240      DO 18 l=1,llm-1
241         z2=-h*log(sig(l+1))
242         write(*,'(i5,3x,4f8.4)') l,-h*log(s(l))/kappa,z2-z1,z2
243     &    ,dsig(l)
244         write(14,'(3x,i5,1f10.4)') l,-h*log(s(l))/kappa
245         z1=z2
24618    CONTINUE
247      write(*,'(i5,3x,3f8.4)') l,-h*log(s(llm))/kappa
248      write(14,'(3x,i5,1f10.4)') l,-h*log(s(llm))/kappa
249
250
251      RETURN
252      END
253
254c---------------------------------------------------------------
255      subroutine adv_vert2(w,q,plev)
256c   Schéma amont pour l'advection verticale de l'eau
257      implicit none
258
259#include "dimensions.h"
260#include "dimphy.h"
261
262      real w(llm+1),wq(llm+1),q(llm),plev(llm+1),dsig(llm+1)
263
264      integer l
265
266      do l=1,llm
267         wq(l)=q(l)*w(l)
268         dsig(l)=(plev(l)-plev(l+1))/plev(1)
269      enddo
270
271      wq(llm+1)=0.
272      do l=1,llm
273         q(l)=(q(l)*dsig(l)+wq(l+1)-wq(l))/(dsig(l)+w(l+1)-w(l))
274      enddo
275      do l=1,llm
276         wq(l)=q(l)*w(l)
277      enddo
278      wq(llm+1)=0.
279      do l=1,llm
280         q(l)=(q(l)*dsig(l)+wq(l+1)-wq(l))/(dsig(l)+w(l+1)-w(l))
281      enddo
282
283      return
284      end
285c---------------------------------------------------------------
286      subroutine adv_vert(w,q)
287c   Schéma amont pour l'advection verticale de l'eau
288      implicit none
289
290#include "dimensions.h"
291#include "dimphy.h"
292#include "comvert1d.h"
293
294      real w(llm+1),wq(llm+1),q(llm)
295
296      integer l
297
298      do l=1,llm
299         wq(l)=q(l)*w(l)
300      enddo
301
302      wq(llm+1)=0.
303      do l=1,llm
304         q(l)=(q(l)*dsig(l)+wq(l+1)-wq(l))/(dsig(l)+w(l+1)-w(l))
305      enddo
306      do l=1,llm
307         wq(l)=q(l)*w(l)
308      enddo
309      wq(llm+1)=0.
310      do l=1,llm
311         q(l)=(q(l)*dsig(l)+wq(l+1)-wq(l))/(dsig(l)+w(l+1)-w(l))
312      enddo
313
314      return
315      end
316      subroutine ug_vg(lm,sig,ug,vg)
317      implicit none
318
319      integer lm
320
321      real u0,uinf,sigu0,alphau
322      real v0,vinf,sigv0,alphav
323
324      real sig(lm),ug(lm),vg(lm)
325      integer l
326
327      data u0,uinf,sigu0/-9.,6.,0.8/
328c     data v0,vinf,sigv0/-1.,4.,0.98/
329      data v0,vinf,sigv0/-1.,7.,0.98/
330
331
332      alphau=-log((uinf-u0)/uinf)/log(sigu0)
333      alphav=-log((vinf-v0)/vinf)/log(sigv0)
334      print*,'Alpha=',alphau,alphav
335
336      do l=1,lm
337         ug(l)=uinf+(u0-uinf)*sig(l)**alphau
338         vg(l)=vinf+(v0-vinf)*sig(l)**alphav
339      enddo
340
341      return
342      end
343      SUBROUTINE abort_gcm(modname, message, ierr)
344     
345      USE IOIPSL
346C
347C Stops the simulation cleanly, closing files and printing various
348C comments
349C
350Input: modname = name of calling program
351C         message = stuff to print
352C         ierr    = severity of situation ( = 0 normal )
353
354      character*20 modname
355      integer ierr
356      character*80 message
357
358      write(*,*) 'in abort_gcm'
359      call histclo
360c     call histclo(2)
361c     call histclo(3)
362c     call histclo(4)
363c     call histclo(5)
364      write(*,*) 'out of histclo'
365      write(*,*) 'Stopping in ', modname
366      write(*,*) 'Reason = ',message
367      if (ierr .eq. 0) then
368        write(*,*) 'Everything is cool'
369      else
370        write(*,*) 'Houston, we have a problem ', ierr
371      endif
372      STOP
373      END
374      FUNCTION q_sat(kelvin, millibar)
375c
376      IMPLICIT none
377c======================================================================
378c Autheur(s): Z.X. Li (LMD/CNRS)
379c Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
380c======================================================================
381c Arguments:
382c kelvin---input-R: temperature en Kelvin
383c millibar--input-R: pression en mb
384c
385c q_sat----output-R: vapeur d'eau saturante en kg/kg
386c======================================================================
387c
388      REAL q_sat, kelvin, millibar
389c
390      REAL r2es
391      PARAMETER (r2es=611.14 *18.0153/28.9644)
392c
393      REAL r3les, r3ies, r3es
394      PARAMETER (R3LES=17.269)
395      PARAMETER (R3IES=21.875)
396c
397      REAL r4les, r4ies, r4es
398      PARAMETER (R4LES=35.86)
399      PARAMETER (R4IES=7.66)
400c
401      REAL rtt
402      PARAMETER (rtt=273.16)
403c
404      REAL retv
405      PARAMETER (retv=28.9644/18.0153 - 1.0)
406c
407      REAL zqsat
408      REAL temp, pres
409C     ------------------------------------------------------------------
410c
411c
412      temp = kelvin
413      pres = millibar * 100.0
414c      write(*,*)'kelvin,millibar=',kelvin,millibar
415c      write(*,*)'temp,pres=',temp,pres
416c
417      IF (temp .LE. rtt) THEN
418         r3es = r3ies
419         r4es = r4ies
420      ELSE
421         r3es = r3les
422         r4es = r4les
423      ENDIF
424c
425      zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
426      zqsat=MIN(0.5,ZQSAT)
427      zqsat=zqsat/(1.-retv  *zqsat)
428c
429      q_sat = zqsat
430c
431      RETURN
432      END
433      subroutine scopy(n,sx,incx,sy,incy)
434c
435      IMPLICIT NONE
436c
437      integer n,incx,incy,ix,iy,i
438      real sx((n-1)*incx+1),sy((n-1)*incy+1)
439c
440      iy=1
441      ix=1
442      do 10 i=1,n
443         sy(iy)=sx(ix)
444         ix=ix+incx
445         iy=iy+incy
44610    continue
447c
448      return
449      end
450      subroutine wrgradsfi(if,nl,field,name,titlevar)
451      implicit none
452
453c   Declarations
454
455#include "gradsdef.h"
456
457c   arguments
458      integer if,nl
459      real field(imx*jmx*lmx)
460      character*10 name,file
461      character*10 titlevar
462
463c   local
464
465      integer im,jm,lm,i,j,l,lnblnk,iv,iii,iji,iif,ijf
466
467      logical writectl
468
469
470      writectl=.false.
471
472c     print*,if,iid(if),jid(if),ifd(if),jfd(if)
473      iii=iid(if)
474      iji=jid(if)
475      iif=ifd(if)
476      ijf=jfd(if)
477      im=iif-iii+1
478      jm=ijf-iji+1
479      lm=lmd(if)
480
481c     print*,'im,jm,lm,name,firsttime(if)'
482c     print*,im,jm,lm,name,firsttime(if)
483
484      if(firsttime(if)) then
485         if(name.eq.var(1,if)) then
486            firsttime(if)=.false.
487            ivar(if)=1
488         print*,'fin de l initialiation de l ecriture du fichier'
489         print*,file
490           print*,'fichier no: ',if
491           print*,'unit ',unit(if)
492           print*,'nvar  ',nvar(if)
493           print*,'vars ',(var(iv,if),iv=1,nvar(if))
494         else
495            ivar(if)=ivar(if)+1
496            nvar(if)=ivar(if)
497            var(ivar(if),if)=name
498            tvar(ivar(if),if)=titlevar(1:lnblnk(titlevar))
499            nld(ivar(if),if)=nl
500            print*,'initialisation ecriture de ',var(ivar(if),if)
501            print*,'if ivar(if) nld ',if,ivar(if),nld(ivar(if),if)
502         endif
503         writectl=.true.
504         itime(if)=1
505      else
506         ivar(if)=mod(ivar(if),nvar(if))+1
507         if (ivar(if).eq.nvar(if)) then
508            writectl=.true.
509            itime(if)=itime(if)+1
510         endif
511
512         if(var(ivar(if),if).ne.name) then
513           print*,'Il faut stoker la meme succession de champs a chaque'
514           print*,'pas de temps'
515           print*,'fichier no: ',if
516           print*,'unit ',unit(if)
517           print*,'nvar  ',nvar(if)
518           print*,'vars ',(var(iv,if),iv=1,nvar(if))
519
520           stop
521         endif
522      endif
523
524c     print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
525c     print*,ivar(if),nvar(if),var(ivar(if),if),writectl
526      do l=1,nl
527         irec(if)=irec(if)+1
528c        print*,'Ecrit rec=',irec(if),iii,iif,iji,ijf,
529c    s (l-1)*imd(if)*jmd(if)+(iji-1)*imd(if)+iii
530c    s ,(l-1)*imd(if)*jmd(if)+(ijf-1)*imd(if)+iif
531         write(unit(if)+1,rec=irec(if))
532     s   ((field((l-1)*imd(if)*jmd(if)+(j-1)*imd(if)+i)
533     s   ,i=iii,iif),j=iji,ijf)
534      enddo
535      if (writectl) then
536
537      file=fichier(if)
538c   WARNING! on reecrase le fichier .ctl a chaque ecriture
539      open(unit(if),file=file(1:lnblnk(file))//'.ctl',
540     &         form='formatted',status='unknown')
541      write(unit(if),'(a5,1x,a40)')
542     &       'DSET ','^'//file(1:lnblnk(file))//'.dat'
543
544      write(unit(if),'(a12)') 'UNDEF 1.0E30'
545      write(unit(if),'(a5,1x,a40)') 'TITLE ',title(if)
546      call formcoord(unit(if),im,xd(iii,if),1.,.false.,'XDEF')
547      call formcoord(unit(if),jm,yd(iji,if),1.,.true.,'YDEF')
548      call formcoord(unit(if),lm,zd(1,if),1.,.false.,'ZDEF')
549      write(unit(if),'(a4,i10,a30)')
550     &       'TDEF ',itime(if),' LINEAR 07AUG1998 30MN '
551      write(unit(if),'(a4,2x,i5)') 'VARS',nvar(if)
552      do iv=1,nvar(if)
553c        print*,'if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)'
554c        print*,if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)
555         write(unit(if),1000) var(iv,if),nld(iv,if)-1/nld(iv,if)
556     &     ,99,tvar(iv,if)
557      enddo
558      write(unit(if),'(a7)') 'ENDVARS'
559c
5601000  format(a5,3x,i4,i3,1x,a39)
561
562      close(unit(if))
563
564      endif ! writectl
565
566      return
567
568      END
569
570      subroutine inigrads(if,im
571     s  ,x,fx,xmin,xmax,jm,y,ymin,ymax,fy,lm,z,fz
572     s  ,dt,file,titlel)
573
574
575      implicit none
576
577      integer if,im,jm,lm,i,j,l,lnblnk
578      real x(im),y(jm),z(lm),fx,fy,fz,dt
579      real xmin,xmax,ymin,ymax
580
581      character file*10,titlel*40
582
583#include "gradsdef.h"
584
585      data unit/24,32,34,36,38,40,42,44,46,48/
586      data nf/0/
587
588      if (if.le.nf) stop'verifier les appels a inigrads'
589
590      print*,'Entree dans inigrads'
591
592      nf=if
593      title(if)=titlel
594      ivar(if)=0
595
596      fichier(if)=file(1:lnblnk(file))
597
598      firsttime(if)=.true.
599      dtime(if)=dt
600
601      iid(if)=1
602      ifd(if)=im
603      imd(if)=im
604      do i=1,im
605         xd(i,if)=x(i)*fx
606         if(xd(i,if).lt.xmin) iid(if)=i+1
607         if(xd(i,if).le.xmax) ifd(if)=i
608      enddo
609      print*,'On stoke du point ',iid(if),'  a ',ifd(if),' en x'
610
611      jid(if)=1
612      jfd(if)=jm
613      jmd(if)=jm
614      do j=1,jm
615         yd(j,if)=y(j)*fy
616         if(yd(j,if).gt.ymax) jid(if)=j+1
617         if(yd(j,if).ge.ymin) jfd(if)=j
618      enddo
619      print*,'On stoke du point ',jid(if),'  a ',jfd(if),' en y'
620
621      print*,'Open de dat'
622      print*,'file=',file
623      print*,'fichier(if)=',fichier(if)
624
625      print*,4*(ifd(if)-iid(if))*(jfd(if)-jid(if))
626      print*,file(1:lnblnk(file))//'.dat'
627
628      OPEN (unit(if)+1,FILE=file(1:lnblnk(file))//'.dat',
629     s   FORM='UNFORMATTED',
630     s   ACCESS='DIRECT'
631     s  ,RECL=4*(ifd(if)-iid(if)+1)*(jfd(if)-jid(if)+1))
632
633      print*,'Open de dat ok'
634
635      lmd(if)=lm
636      do l=1,lm
637         zd(l,if)=z(l)*fz
638      enddo
639
640      irec(if)=0
641
642      print*,if,imd(if),jmd(if),lmd(if)
643      print*,'if,imd(if),jmd(if),lmd(if)'
644
645      return
646      end
647      SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
648      IMPLICIT NONE
649c=======================================================================
650c   passage d'un champ de la grille scalaire a la grille physique
651c=======================================================================
652
653c-----------------------------------------------------------------------
654c   declarations:
655c   -------------
656
657      INTEGER im,jm,ngrid,nfield
658      REAL pdyn(im,jm,nfield)
659      REAL pfi(ngrid,nfield)
660
661      INTEGER j,ifield,ig
662
663c-----------------------------------------------------------------------
664c   calcul:
665c   -------
666
667      IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)
668     s    STOP 'probleme de dim'
669c   traitement des poles
670      CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
671      CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
672
673c   traitement des point normaux
674      DO ifield=1,nfield
675         DO j=2,jm-1
676            ig=2+(j-2)*(im-1)
677            CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
678         ENDDO
679      ENDDO
680
681      RETURN
682      END
683
684      subroutine writelim
685     s   (phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice)
686c
687#include "dimensions.h"
688#include "dimphy.h"
689#include "netcdf.inc"
690
691      REAL phy_nat(klon,360)
692      REAL phy_alb(klon,360)
693      REAL phy_sst(klon,360)
694      REAL phy_bil(klon,360)
695      REAL phy_rug(klon,360)
696      REAL phy_ice(klon,360)
697
698      INTEGER ierr
699      INTEGER dimfirst(3)
700      INTEGER dimlast(3)
701c
702      INTEGER nid, ndim, ntim
703      INTEGER dims(2), debut(2), epais(2)
704      INTEGER id_tim
705      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
706
707      PRINT*, 'Ecriture du fichier limit'
708c
709      ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
710c
711      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
712     .                       "Fichier conditions aux limites")
713      ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
714      ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
715c
716      dims(1) = ndim
717      dims(2) = ntim
718c
719ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
720      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
721      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
722     .                        "Jour dans l annee")
723ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
724      ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
725      ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
726     .                        "Nature du sol (0,1,2,3)")
727ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
728      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
729      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
730     .                        "Temperature superficielle de la mer")
731ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
732      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
733      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
734     .                        "Reference flux de chaleur au sol")
735ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
736      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
737      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
738     .                        "Albedo a la surface")
739ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
740      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
741      ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
742     .                        "Rugosite")
743c
744      ierr = NF_ENDDEF(nid)
745c
746      DO k = 1, 360
747c
748      debut(1) = 1
749      debut(2) = k
750      epais(1) = klon
751      epais(2) = 1
752c
753#ifdef NC_DOUBLE
754      ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
755      ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
756      ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
757      ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
758      ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
759      ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
760#else
761      ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
762      ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
763      ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
764      ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
765      ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
766      ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
767#endif
768c
769      ENDDO
770c
771      ierr = NF_CLOSE(nid)
772c
773      return
774      end
775
776      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
777
778c    Auteur :  P. Le Van .
779c
780      IMPLICIT NONE
781
782#include "dimensions.h"
783#include "paramet.h"
784c
785c=======================================================================
786c
787c
788c    s = sigma ** kappa   :  coordonnee  verticale
789c    dsig(l)            : epaisseur de la couche l ds la coord.  s
790c    sig(l)             : sigma a l'interface des couches l et l-1
791c    ds(l)              : distance entre les couches l et l-1 en coord.s
792c
793c=======================================================================
794c
795      REAL pa,preff
796      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
797      REAL presnivs(llm)
798c
799c   declarations:
800c   -------------
801c
802      REAL sig(llm+1),dsig(llm)
803c
804      INTEGER l
805      REAL snorm
806      REAL alpha,beta,gama,delta,deltaz,h
807      INTEGER np,ierr
808      REAL pi,x
809
810c-----------------------------------------------------------------------
811c
812      pi=2.*ASIN(1.)
813
814      OPEN(99,file='sigma.def',status='old',form='formatted',
815     s   iostat=ierr)
816
817c-----------------------------------------------------------------------
818c   cas 1 on lit les options dans sigma.def:
819c   ----------------------------------------
820
821      IF (ierr.eq.0) THEN
822
823      print*,'WARNING!!! on lit les options dans sigma.def'
824      READ(99,*) deltaz
825      READ(99,*) h
826      READ(99,*) beta
827      READ(99,*) gama
828      READ(99,*) delta
829      READ(99,*) np
830      CLOSE(99)
831      alpha=deltaz/(llm*h)
832c
833
834       DO 1  l = 1, llm
835       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*
836     $          ( (tanh(gama*l)/tanh(gama*llm))**np +
837     $            (1.-l/FLOAT(llm))*delta )
838   1   CONTINUE
839
840       sig(1)=1.
841       DO 101 l=1,llm-1
842          sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
843101    CONTINUE
844       sig(llm+1)=0.
845
846       DO 2  l = 1, llm
847       dsig(l) = sig(l)-sig(l+1)
848   2   CONTINUE
849c
850
851      ELSE
852c-----------------------------------------------------------------------
853c   cas 2 ancienne discretisation (LMD5...):
854c   ----------------------------------------
855
856      PRINT*,'WARNING!!! Ancienne discretisation verticale'
857
858      h=7.
859      snorm  = 0.
860      DO l = 1, llm
861         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
862         dsig(l) = 1.0 + 7.0 * SIN(x)**2
863         snorm = snorm + dsig(l)
864      ENDDO
865      snorm = 1./snorm
866      DO l = 1, llm
867         dsig(l) = dsig(l)*snorm
868      ENDDO
869      sig(llm+1) = 0.
870      DO l = llm, 1, -1
871         sig(l) = sig(l+1) + dsig(l)
872      ENDDO
873
874      ENDIF
875
876
877      DO l=1,llm
878        nivsigs(l) = FLOAT(l)
879      ENDDO
880
881      DO l=1,llmp1
882        nivsig(l)= FLOAT(l)
883      ENDDO
884
885c
886c    ....  Calculs  de ap(l) et de bp(l)  ....
887c    .........................................
888c
889c
890c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
891c
892
893      bp(llmp1) =   0.
894
895      DO l = 1, llm
896cc
897ccc    ap(l) = 0.
898ccc    bp(l) = sig(l)
899
900      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
901      ap(l) = pa * ( sig(l) - bp(l) )
902c
903      ENDDO
904      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
905
906      PRINT *,' BP '
907      PRINT *,  bp
908      PRINT *,' AP '
909      PRINT *,  ap
910
911      DO l = 1, llm
912       dpres(l) = bp(l) - bp(l+1)
913       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
914      ENDDO
915
916      PRINT *,' PRESNIVS '
917      PRINT *,presnivs
918
919      RETURN
920      END
Note: See TracBrowser for help on using the repository browser.