source: LMDZ.3.3/trunk/libf/phylmd/1Dconv.h @ 5415

Last change on this file since 5415 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: 54.3 KB
Line 
1        SUBROUTINE get_uvd(itap,dtime,tsol,qsol,file_fordat
2     s                    ,ht,hq,hw)
3
4        implicit none
5 
6cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
7c cette routine permet d'obtenir u_convg,v_convg,ht,hq et ainsi de
8c pouvoir calculer la convergence et le cisaillement dans la physiq
9ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
10
11#include "YOMCST.h"
12
13      INTEGER klev
14      REAL play(100)  !pression en Pa au milieu de chaque couche GCM
15      INTEGER JM(100) !pression en Pa au milieu de chaque couche GCM
16      REAL coef1(100) !coefficient d'interpolation
17      REAL coef2(100) !coefficient d'interpolation
18
19      INTEGER nblvlm !nombre de niveau de pression du mesoNH
20      REAL playm(100)  !pression en Pa au milieu de chaque couche Meso-NH
21      REAL hplaym(100) !pression en hPa milieux des couches Meso-NH
22
23      integer i,j,k,ii,ll,in
24      REAL tsol,qsol
25
26      CHARACTER*80 file_forctl,file_fordat,file_start
27
28      COMMON/physiq1/klev,play,JM,coef1,coef2
29      COMMON/physiq2/nblvlm,playm,hplaym
30
31c======================================================================
32c methode: on va chercher les donnees du mesoNH de meteo france, on y
33c          a acces a tout pas detemps grace a la routine rdgrads qui
34c          est une boucle lisant dans ces fichiers.
35c          Puis on interpole ces donnes sur les 11 niveaux du gcm et
36c          et sur les pas de temps de ce meme gcm
37c======================================================================
38c input:
39c       pasmax     :nombre de pas de temps maximum du mesoNH
40c       dt         :pas de temps du meso_NH (en secondes)
41c----------------------------------------------------------------------
42      integer pasmax,dt
43      save pasmax,dt
44c----------------------------------------------------------------------
45c arguments:
46c           itap   :compteur de la physique(le nombre de ces pas est
47c                   fixe dans la subroutine calcul_ini_gcm de interpo
48c                   -lation
49c           dtime  :pas detemps du gcm (en secondes)
50c           ht     :convergence horizontale de temperature(K/s)
51c           hq     :    "         "       d'humidite (kg/kg/s)
52c           hw     :vitesse verticale moyenne (m/s**2)
53c----------------------------------------------------------------------
54        integer itap
55        real dtime
56        real ht(100)
57        real hq(100)
58        real hw(100)
59c----------------------------------------------------------------------
60c Variables internes de get_uvd (note : l'interpolation temporelle
61c est faite entre les pas de temps before et after, sur les variables
62c definies sur la grille du SCM)
63c     time0     :date initiale en secondes
64c     time      :temps associe a chaque pas
65c     pas       :numero du pas du meso_NH
66c     pasprev   :numero du pas precedent
67c     htaft     :advection horizontale de temp. au pas de temps after
68c     hqaft     :    "         "      d'humidite        "
69c     hwaft     :vitesse verticalle moyenne  au pas de temps after
70c     htbef     :idem htaft, mais pour le pas de temps before
71c     hqbef     :voir hqaft
72c     hwbef     :voir hwaft
73c----------------------------------------------------------------------
74        integer time0,pas,pasprev
75        save time0,pas,pasprev
76        real time
77        real htaft(100),hqaft(100),hwaft(100)
78        save htaft,hqaft,hwaft
79        real htbef(100),hqbef(100),hwbef(100)
80        save htbef,hqbef,hwbef
81        integer timeaft,timebef
82        save timeaft,timebef
83        integer temps
84        character*4 string
85c----------------------------------------------------------------------
86c variables arguments de la subroutine rdgrads
87c---------------------------------------------------------------------
88        integer icompt      !compteur de rdgrads
89        real z(100)         ! altitude (grille Meso)
90        real ht_mes(100)    !convergence horizontale de temperature
91                            !-(grille Meso)
92        real hq_mes(100)    !convergence horizontale d'humidite
93                            !(grille Meso)
94        real hw_mes(100)    !vitesse verticale moyenne
95                            !(grille Meso)
96c
97c---------------------------------------------------------------------
98c variable argument de la subroutine copie
99c---------------------------------------------------------------------
100c SB        real pplay(100)    !pression en milieu de couche du gcm
101c SB                            !argument de la physique
102c---------------------------------------------------------------------
103c variables destinees a la lecture du pas de temps du fichier de donnees
104c---------------------------------------------------------------------
105       character*80 aaa,atemps,spaces,apasmax
106       integer nch,imn,ipa
107c---------------------------------------------------------------------
108c  procedures appelees
109        external rdgrads    !lire en iterant dans forcing.dat
110c---------------------------------------------------------------------
111               print*,'le pas itap est:',itap
112c*** on determine le pas du meso_NH correspondant au nouvel itap ***
113c*** pour aller chercher les champs dans rdgrads                 ***
114        time=time0+itap*dtime
115        temps=int(time/dt+1)
116        pas=min(temps,pasmax)
117             print*,'le pas Meso est:',pas
118c
119c
120c===================================================================
121c
122c*** on remplit les champs before avec les champs after du pas   ***
123c*** precedent en format gcm                                     ***
124        if(pas.gt.pasprev)then
125               do i=1,klev
126                  htbef(i)=htaft(i)
127                  hqbef(i)=hqaft(i)
128                  hwbef(i)=hwaft(i)
129               enddo
130               timebef=pasprev*dt
131               timeaft=timebef+dt
132               icompt=(pas-1)*(nblvlm*4)
133                       print*,'le pas pas est:',pas
134c*** on va chercher les nouveaux champs after dans toga.dat     ***
135c*** champs en format meso_NH                                   ***
136c          open(99,FILE='forcing.dat',FORM='UNFORMATTED',
137
138          write(*,'(a)') 'OPEN dans get_uvd de '//file_fordat
139          open(99,FILE=file_fordat,FORM='UNFORMATTED',
140     .             ACCESS='DIRECT',RECL=4)
141          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes)
142          do i = 1,nblvlm
143            ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
144          enddo
145c
146               print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
147               print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
148               print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
149c*** on interpole les champs meso_NH sur les niveaux de pression***
150c*** gcm . on obtient le nouveau champ after                    ***
151            do k=1,klev
152             if (JM(k) .eq. 0) then
153         htaft(k)=coef1(k)*tsol+coef2(k)*ht_mes(jm(k)+1)
154         hqaft(k)=coef1(k)*qsol+coef2(k)*hq_mes(jm(k)+1)
155         hwaft(k)=              coef2(k)*hw_mes(jm(k)+1)
156             else
157         htaft(k)=coef1(k)*ht_mes(jm(k))+coef2(k)*ht_mes(jm(k)+1)
158         hqaft(k)=coef1(k)*hq_mes(jm(k))+coef2(k)*hq_mes(jm(k)+1)
159         hwaft(k)=coef1(k)*hw_mes(jm(k))+coef2(k)*hw_mes(jm(k)+1)
160             endif
161            enddo
162          pasprev=pas
163         else
164                      print*,'timebef est:',timebef
165         endif      !fin du bloc relatif au passage au pas
166                    !de temps (meso) suivant
167c*** si on atteint le pas max des donnees experimentales ,on     ***
168c*** on conserve les derniers champs calcules                    ***
169      if(pas.ge.pasmax)then
170          do ll=1,klev
171               ht(ll)=htaft(ll)
172               hq(ll)=hqaft(ll)
173               hw(ll)=hwaft(ll)
174          enddo
175      else
176c*** on interpole sur les pas de temps de 10mn du gcm a partir   ***
177c** des pas de temps de 1h du meso_NH                            ***
178         do j=1,klev
179         ht(j)=((timeaft-time)*htbef(j)+(time-timebef)*htaft(j))/dt
180         hq(j)=((timeaft-time)*hqbef(j)+(time-timebef)*hqaft(j))/dt
181         hw(j)=((timeaft-time)*hwbef(j)+(time-timebef)*hwaft(j))/dt
182         enddo
183       endif
184c
185c-------------------------------------------------------------------
186c
187         return
188c
189c-----------------------------------------------------------------------
190c on sort les champs de "convergence" pour l'instant initial 'in'
191c ceci se passe au pas temps itap=0 de la physique
192c-----------------------------------------------------------------------
193        entry get_uvd2(itap,file_forctl,file_fordat,file_start
194     s                ,ht,hq,hw)
195             print*,'le pas itap est:',itap
196c
197c===================================================================
198c
199      write(*,*) '   '
200      write(*,*) 'FICHIERS A LIRE DANS GET_UVD2:   '
201      write(*,'(a)') 'fichier forcing.ctl: '//file_forctl
202      write(*,'(a)') 'fichier forcing.dat: '//file_fordat
203      write(*,'(a)') 'fichier start18.data: '//file_start
204      write(*,*) '   '
205
206c!! en attendant de pouvoir compiler les fns CERN, en prescrit
207c!! les variables imn et pasmax a la main...
208c!!
209       write(*,'(a)') 'OPEN '//file_forctl
210       open(97,FILE=file_forctl,FORM='FORMATTED')
211c
212c------------------
213      do i=1,1000
214      read(97,1000,end=999) string
215 1000 format (a4)
216      if (string .eq. 'TDEF') go to 50
217      enddo
218 50   backspace(97)
219c-------------------------------------------------------------------
220c   *** on lit le pas de temps dans le fichier de donnees ***
221c   *** "forcing.ctl" et pasmax                           ***
222c-------------------------------------------------------------------
223      read(97,2000) aaa
2242000  format (a80)
225         print*,'aaa est',aaa
226      aaa=spaces(aaa,1)
227         print*,'aaa',aaa
228      call getsch(aaa,' ',' ',5,atemps,nch)
229         print*,'atemps est',atemps
230        atemps=atemps(1:nch-2)
231         print*,'atemps',atemps
232        read(atemps,*) imn
233        dt=imn*60
234         print*,'le pas de temps dt',dt
235      call getsch(aaa,' ',' ',2,apasmax,nch)
236        apasmax=apasmax(1:nch)
237        read(apasmax,*) ipa
238        pasmax=ipa
239         print*,'pasmax est',pasmax
240      CLOSE(97)
241
242c CASE_E:
243c!!         imn = 60
244c!!         ipa = 8
245c TOGA:
246c!!         imn = 360
247c!!         ipa = 480
248
249         dt=imn*60
250         pasmax=ipa
251         print*,'le pas de temps dt',dt
252         print*,'pasmax est',pasmax
253
254
255c------------------------------------------------------------------
256c *** onlit le pas de temps initial de la simulation dans ***
257c *** "start.data"                                        ***
258c------------------------------------------------------------------
259c      open(98,file='start18.data',form='formatted')
260      write(*,'(a)') 'OPEN '//file_start
261      open(98,FILE=file_start,FORM='FORMATTED')
262          read(98,*)in
263                  pasprev=in
264                    print*,'le pas in ini est:',pasprev
265C
266Cjyg     Correction de la date du demarrage.
267CC                  time0=dt*pasprev
268                  time0=dt*(pasprev-1)
269C
270          close(98)
271c
272c      open(99,FILE='forcing.dat',FORM='UNFORMATTED',
273      write(*,'(a)') 'OPEN '//file_fordat
274      open(99,FILE=file_fordat,FORM='UNFORMATTED',
275     .          ACCESS='DIRECT',RECL=4)
276                  icompt=(in-1)*(nblvlm*4)
277        call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes)
278          do i = 1,nblvlm
279            ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
280          enddo
281c
282               print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
283               print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
284               print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
285c----------------------------------------------------------------------
286c on a obtenu des champs initiaux sur les niveaux du meso_NH
287c on interpole sur les niveaux du gcm(niveau pression bien sur!)
288c-----------------------------------------------------------------------
289         do ii=1,klev
290         htaft(ii)=coef1(ii)*ht_mes(JM(ii))+coef2(ii)*ht_mes(JM(ii)+1)
291         hqaft(ii)=coef1(ii)*hq_mes(JM(ii))+coef2(ii)*hq_mes(JM(ii)+1)
292         hwaft(ii)=coef1(ii)*hw_mes(JM(ii))+coef2(ii)*hw_mes(JM(ii)+1)
293             enddo
294c valeurs initiales des champs de convergence
295          do k=1,klev
296             ht(k)=htaft(k)
297             hq(k)=hqaft(k)
298             hw(k)=hwaft(k)
299          enddo
300        close(99)
301        close(98)
302c
303c-------------------------------------------------------------------
304c
305 100    return
306c
307999     continue
308        stop 'erreur lecture, file forcing.ctl'
309        end
310 
311 
312      SUBROUTINE cool_pool(istep
313     e                    ,n_cooling,dt_cooling,dq_cooling
314     s                    ,dt_cool,dq_cool)
315       implicit none
316C***************************************************************
317C*                                                             *
318C* COOL_POOL                                                   *
319C*                                                             *
320C*                                                             *
321C* written by   : Gilles Foret RAMSES, 15/09/97, 22.00.2       *
322C* modified by :  Sandrine Bony 10/09/98                       *
323C***************************************************************
324c Arguments
325c =========
326c Input
327c -----
328c   istep : numero du pas de temps
329c   n_cooling: nbre de pas de temps ou la pertubation nominale
330c              est appliquee (ensuite, la pertubation decroit
331c              exponentiellement).
332c   dt_cooling : pertubation nominale en temperature
333c   dq_cooling : pertubation nominale en humidite
334c Output
335c ------
336c    dt_cool : pertubation en temperature
337c    dq_cool : pertubation en humidite
338c
339c Variables internes
340c ==================
341c    scale :  facteur applique a la pertubation nominale
342c
343#include "dimensions.h"
344#include "dimphy.h"
345c
346      integer n_cooling,k,istep
347      real dt_cooling(klev),dq_cooling(klev),scale
348      real dt_cool(klev),dq_cool(klev)
349c
350      if (istep .le. n_cooling ) then
351        scale = 1.
352      else
353        scale = 4**(min(15,istep-n_cooling))
354      endif
355c
356        do k = 1,klev
357          dt_cool(k) = dt_cooling(k)/scale
358          dq_cool(k) = dq_cooling(k)/scale
359        enddo
360c
361      return
362      end
363      SUBROUTINE advect_tvl(dtime,t,q,vu_f,vv_f,t_f,q_f
364     :                     ,d_t_adv,d_q_adv)
365      implicit none
366
367#include "dimensions.h"
368#include "dimphy.h"
369
370      integer k
371      real dtime, fact, du, dv, cx, cy, alx, aly
372      real t(klev), q(klev,3)
373     :   , vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
374
375      real d_t_adv(klev), d_q_adv(klev,3)
376
377c Velocity of moving cell
378      data cx,cy /12., -2./
379
380c Dimensions of moving cell
381      data alx,aly /100 000.,150 000./
382
383      do k = 1, klev
384            du = abs(vu_f(k)-cx)/alx
385            dv = abs(vv_f(k)-cy)/aly
386            fact = dtime *(du+dv-du*dv*dtime)
387            d_t_adv(k) = fact * (t_f(k)-t(k))
388            d_q_adv(k,1) = fact * (q_f(k,1)-q(k,1))
389            d_q_adv(k,2) = fact * (q_f(k,2)-q(k,2))
390            d_q_adv(k,3) = fact * (q_f(k,3)-q(k,3))
391      enddo
392
393      return
394      end
395      SUBROUTINE copie(klevgcm,playgcm,psolgcm,file_forctl)
396      implicit none
397
398ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
399c cette routine remplit les COMMON physiq1 et physiq2.h
400cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
401
402      INTEGER JM
403      INTEGER klev !nombre de niveau de pression du GCM
404      INTEGER nblvlm !nombre de niveau de pression du mesoNH
405
406      REAL playm(100) !pression en Pa au milieu de chaque couche Meso-NH
407      REAL hplaym(100)!pression en hecto-Pa des milieux de couche Meso-NH
408      REAL play(100)  !pression en Pa au milieu de chaque couche GCM
409      REAL coef1(100)   !coefficient d'interpolation
410      REAL coef2(100)   !coefficient d'interpolation
411
412      COMMON/physiq1/klev,play,JM,coef1,coef2
413      COMMON/physiq2/nblvlm,playm,hplaym
414
415      integer i,k,klevgcm
416      real playgcm(klevgcm) ! pression en milieu de couche du gcm
417      real psolgcm
418      character*80 file_forctl
419
420      klev = klevgcm
421
422c---------------------------------------------------------------------
423c pression au milieu des couches du gcm dans la physiq
424c (SB: remplace le call conv_lipress_gcm(playgcm) )
425c---------------------------------------------------------------------
426
427       do k = 1, klev
428        play(k) = playgcm(k)
429        print*,'la pression gcm est:',play(k)
430       enddo
431
432c----------------------------------------------------------------------
433c lecture du descripteur des donnees Meso-NH (forcing.ctl):
434-> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
435c (on remplit le COMMON physiq2)
436c----------------------------------------------------------------------
437
438      call mesolupbis(file_forctl)
439
440      print*,'la valeur de nblvlm est:',nblvlm
441
442c----------------------------------------------------------------------
443c etude de la correspondance entre les niveaux meso.NH et GCM;
444c calcul des coefficients d'interpolation coef1 et coef2
445c (on remplit le COMMON physiq1)
446c----------------------------------------------------------------------
447
448      call corresbis(psolgcm)
449
450c---------------------------------------------------------
451c TEST sur le remplissage de physiq1 et physiq2:
452c---------------------------------------------------------
453 
454      write(*,*) ' '
455      write(*,*) 'TESTS physiq1 et physiq2 dans copie.F '
456      write(*,*) '--------------------------------------'
457      write(*,*) 'GCM: nb niveaux:',klev,' et pression, coeffs:'
458      do k = 1, klev
459      write(*,*) play(k), coef1(k), coef2(k)
460      enddo
461      write(*,*) 'MESO-NH: nb niveaux:',nblvlm,' et pression:'
462      do k = 1, nblvlm
463      write(*,*) playm(k), hplaym(k)
464      enddo
465      write(*,*) ' '
466
467      end
468      SUBROUTINE writeg1d(ngrid,nx,x,nom,titre)
469      IMPLICIT NONE
470c.......................................................................
471c
472c  ecriture de x pour GRADS-1D
473c
474in :
475c         * ngrid      ---> pour controler que l'on est bien en 1D
476c         * nx         ---> taille du vecteur a stocker
477c                             "1" pour une variable de surface
478c                             "nlayer" pour une variable de centre de couche
479c                             "nlayer+1" pour une variable d'interface
480c         * x          ---> variable a stocker
481c         * nom        ---> nom "pour grads"
482c         * titre      ---> titre "pour grads"
483c
484c.......................................................................
485c
486#include "comg1d.h"
487c
488c.......................................................................
489c  declaration des arguments
490c
491      INTEGER ngrid,nx
492      REAL x(nx)
493      CHARACTER*(*) nom
494      CHARACTER*(*) titre
495c
496c  declaration des arguments
497c....................................................................... 
498c  declaration des variables locales
499c
500      INTEGER ilayer,ivar
501      LOGICAL test
502c
503c  declaration des variables locales
504c.......................................................................
505c  contole 1D
506c
507c     print*,'ngrid=',ngrid
508      IF (ngrid.NE.1) return
509c
510c  contole 1D
511c.......................................................................
512c  ouverture du fichier au premier appel
513c
514      IF (g1d_premier) THEN
515        OPEN (g1d_unitfich,FILE=g1d_nomfich
516     &       ,FORM='unformatted',ACCESS='direct',RECL=4)
517        g1d_irec=0
518        g1d_nvar=0
519        g1d_premier=.false.
520      ENDIF
521c
522c  ouverture du fichier au premier appel
523c.......................................................................
524c  pour l'ecriture du fichier ctl
525c
526      test=.true.
527      DO ivar=1,g1d_nvar
528        IF (nom.EQ.g1d_nomvar(ivar)) test=.false.
529      ENDDO
530      IF (test) THEN
531        g1d_nvar=g1d_nvar+1
532        g1d_nomvar(g1d_nvar)=nom
533        g1d_titrevar(g1d_nvar)=titre
534        IF (nx.EQ.1) THEN
535           g1d_dimvar(g1d_nvar)=0
536        ELSEIF (nx.EQ.g1d_nlayer) THEN
537           g1d_dimvar(g1d_nvar)=g1d_nlayer
538        ELSEIF (nx.EQ.g1d_nlayer+1) THEN
539           g1d_dimvar(g1d_nvar)=g1d_nlayer
540        ELSE
541           PRINT *,'._. probleme de dimension dans GRADS-1D ._.'
542        ENDIF
543      ENDIF
544c
545c  pour l'ecriture du fichier ctl
546c.......................................................................
547c  ecriture
548c
549      IF (nx.EQ.1) THEN
550        g1d_irec=g1d_irec+1
551        WRITE(g1d_unitfich,REC=g1d_irec) x(1)
552      ELSE
553        DO ilayer=1,g1d_nlayer
554          g1d_irec=g1d_irec+1
555          WRITE(g1d_unitfich,REC=g1d_irec) x(ilayer)
556        ENDDO
557      ENDIF
558c
559c  ecriture
560c.......................................................................
561c
56210001 CONTINUE
563c
564c.......................................................................
565c
566      RETURN
567      END
568
569
570
571
572
573
574c SB      SUBROUTINE endg1d(ngrid,nlayer,zlayer,ndt)
575      SUBROUTINE endg1d(ngrid,nlayer,player,ndt,dt)
576      IMPLICIT NONE
577c.......................................................................
578c
579c  ecriture du fichier de controle pour GRADS-1D
580c
581in :
582c         * ngrid      ---> pour controler que l'on est bien en 1D
583c         * nlayer     ---> nombre de couches
584c         * zlayer     ---> altitude au centre de chaque couche (km)
585c         * player     ---> pression au centre de chaque couche (hPa)
586c         * ndt        ---> nombre de pas de temps
587c         * dt         ---> valeur du pas de temps (s)
588c
589c.......................................................................
590c
591#include "comg1d.h"
592c
593c.......................................................................
594c  declaration des arguments
595c
596      INTEGER ngrid,nlayer
597c SB      REAL zlayer(nlayer)
598      REAL player(nlayer)
599      INTEGER ndt
600      REAL dt,dtm
601c
602c  declaration des arguments
603c....................................................................... 
604c  declaration des variables locales
605c
606      INTEGER ivar,ilayer
607c
608c  declaration des variables locales
609c.......................................................................
610c  contole 1D
611c
612      IF (ngrid.NE.1) GOTO 10001
613c
614c  contole 1D
615c.......................................................................
616c
617      IF (nlayer.ne.g1d_nlayer) 
618     &   PRINT *,'._. probleme de dimension dans GRADS-1D ._.'
619c
620c.......................................................................
621c
622      CLOSE (g1d_unitfich)
623c
624c.......................................................................
625c
626      dtm = dt/60.
627
628      OPEN (g1d_unitctl,FILE=g1d_nomctl,FORM='formatted'
629     s     ,status='new')
630      WRITE (g1d_unitctl,'(a4,2x,a20)') 'DSET',g1d_nomfich
631      WRITE (g1d_unitctl,'(a5,2x,a20)') 'UNDEF ','1.E+30'
632      WRITE (g1d_unitctl,'(a11)') 'FORMAT YREV'
633      WRITE (g1d_unitctl,'(a5,2x,a30)') 'TITLE ','champs 1D'
634      WRITE (g1d_unitctl,'(a5,i4,a20)') 'XDEF ',1,' LINEAR 0 1'
635      WRITE (g1d_unitctl,'(a5,i4,a20)') 'YDEF ',1,' LINEAR 0 1'
636      WRITE (g1d_unitctl,'(a5,i4,a20)') 'ZDEF ',g1d_nlayer,' LEVELS'
637      WRITE (g1d_unitctl,'(5(1x,f13.5))')
638c SB     &      (zlayer(ilayer),ilayer=1,g1d_nlayer)
639     &      (player(ilayer)/100.,ilayer=1,g1d_nlayer)
640c SB      WRITE (g1d_unitctl,'(a4,2x,i10,a25)')
641c SB     &      'TDEF ',ndt,' LINEAR 02JAN1987 1HR '
642      WRITE (g1d_unitctl,'(a4,2x,i10,a20,i3,a3)')
643     &      'TDEF ',ndt,' LINEAR 02JAN1987 ',INT(dtm),'MN '
644      WRITE (g1d_unitctl,'(a5,i5)') 'VARS ',g1d_nvar
645      DO ivar=1,g1d_nvar
646      WRITE (g1d_unitctl,'(a5,3x,i4,i3,1x,a39)') 
647     &       g1d_nomvar(ivar),g1d_dimvar(ivar),99,g1d_titrevar(ivar)
648      ENDDO
649      WRITE (g1d_unitctl,'(a7)') 'ENDVARS'
650      CLOSE (g1d_unitctl)
651c
652c.......................................................................
653c
65410001 CONTINUE
655c
656c.......................................................................
657c
658      RETURN
659      END
660      SUBROUTINE mesolupbis(file_forctl)
661      implicit none
662c
663cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
664c
665c Lecture descripteur des donnees MESO-NH (forcing.ctl):
666c -------------------------------------------------------
667c
668c     Cette subroutine lit dans le fichier de controle "essai.ctl"
669c     et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
670c     des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
671cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
672c
673      INTEGER nblvlm !nombre de niveau de pression du mesoNH
674      REAL playm(100)  !pression en Pa milieu de chaque couche Meso-NH
675      REAL hplaym(100) !pression en hPa des milieux de couche Meso-NH
676      COMMON/physiq2/nblvlm,playm,hplaym
677
678      INTEGER i,lu,k,mlz,mlzh,j
679 
680      character*80 file_forctl
681
682      character*4 a
683      character*80 aaa,anblvl,spaces
684      integer nch
685
686      lu=9
687c      open (lu,file='forcing.ctl')
688      open(lu,file=file_forctl,form='formatted')
689c
690      do i=1,1000
691      read(lu,1000,end=999) a
692      if (a .eq. 'ZDEF') go to 100
693      enddo
694c
695 100  backspace(lu)
696      print*,'  DESCRIPTION DES 2 MODELES : '
697      print*,' '
698c
699      read(lu,2000) aaa
700 2000  format (a80)
701       aaa=spaces(aaa,1)
702       call getsch(aaa,' ',' ',2,anblvl,nch)
703         read(anblvl,*) nblvlm
704
705c      write(*,*) 'ATTENTION! dans mesolupbis on rentre
706c     : nblvlm a la main car pas de bibliotheque CERN..:'
707c CASE_e:
708c!          nblvlm = 43
709c TOGA:
710c!!          nblvlm = 40
711c
712      print*,'nbre de niveaux de pression Meso-NH :',nblvlm
713      print*,' '
714      print*,'pression en Pa de chaque couche du meso-NH :'
715c
716      read(lu,*) (playm(mlz),mlz=1,nblvlm)
717c      Si la pression est en HPa, la multiplier par 100
718      if (playm(1) .lt. 10000.) then
719        do mlz = 1,nblvlm
720         playm(mlz) = playm(mlz)*100.
721        enddo
722      endif
723      print*,(playm(mlz),mlz=1,nblvlm)
724c
725 1000 format (a4)
726 1001 format(5x,i2)
727c
728      print*,' '
729      do mlzh=1,nblvlm
730      hplaym(mlzh)=playm(mlzh)/100.
731      enddo
732c
733      print*,'pression en hPa de chaque couche du meso-NH: '
734      print*,(hplaym(mlzh),mlzh=1,nblvlm)
735c
736      close (lu)
737      return
738c
739 999  stop 'erreur lecture des niveaux pression des donnees'
740      end
741      SUBROUTINE GETSCH(STR,DEL,TRM,NTH,SST,NCH)
742C***************************************************************
743C*                                                             *
744C*                                                             *
745C* GETSCH                                                      *
746C*                                                             *
747C*                                                             *
748C* modified by :                                               *
749C***************************************************************
750C*   Return in SST the character string found between the NTH-1 and NTH
751C*   occurence of the delimiter 'DEL' but before the terminator 'TRM' in
752C*   the input string 'STR'. If TRM=DEL then STR is considered unlimited.
753C*   NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
754C*   NTH is greater than the number of delimiters in STR.
755      IMPLICIT INTEGER (A-Z)
756      CHARACTER STR*(*),DEL*1,TRM*1,SST*(*)
757      NCH=-1
758      SST=' '
759      IF(NTH.GT.0) THEN
760        IF(TRM.EQ.DEL) THEN
761          LENGTH=LEN(STR)
762        ELSE
763          LENGTH=INDEX(STR,TRM)-1
764          IF(LENGTH.LT.0) LENGTH=LEN(STR)
765        ENDIF
766C*     Find beginning and end of the NTH DEL-limited substring in STR
767        END=-1
768        DO 1,N=1,NTH
769        IF(END.EQ.LENGTH) RETURN
770        BEG=END+2
771        END=BEG+INDEX(STR(BEG:LENGTH),DEL)-2
772        IF(END.EQ.BEG-2) END=LENGTH
773C*        PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
774    1   CONTINUE
775        NCH=END-BEG+1
776        IF(NCH.GT.0) SST=STR(BEG:END)
777      ENDIF
778      END
779      SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw)
780      IMPLICIT none
781
782      INTEGER itape,icount,icomp, nl
783      real z(nl),ht(nl),hq(nl),hw(nl)
784c
785      INTEGER i, k
786c
787      icomp = icount
788c
789c
790         do k=1,nl
791            icomp=icomp+1
792            read(itape,rec=icomp)z(k)
793         enddo
794         do k=1,nl
795            icomp=icomp+1
796            read(itape,rec=icomp)hT(k)
797         enddo
798         do k=1,nl
799            icomp=icomp+1
800            read(itape,rec=icomp)hQ(k)
801         enddo
802         do k=1,nl
803            icomp=icomp+1
804            read(itape,rec=icomp)hw(k)
805         enddo
806c
807c
808      RETURN
809      END
810      SUBROUTINE corresbis(psol)
811      implicit none
812
813ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
814c Cette subroutine calcule et affiche les valeurs des coefficients
815c d'interpolation qui serviront dans la formule d'interpolation elle-
816c meme.
817cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
818
819      INTEGER klev    !nombre de niveau de pression du GCM
820      REAL play(100)  !pression en Pa au milieu de chaque couche GCM
821      INTEGER JM(100)
822      REAL coef1(100) !coefficient d'interpolation
823      REAL coef2(100) !coefficient d'interpolation
824
825      INTEGER nblvlm !nombre de niveau de pression du mesoNH
826      REAL playm(100) !pression en Pa milieu de chaque couche Meso-NH
827      REAL hplaym(100)!pression en hPa des milieux de couche Meso-NH
828
829      COMMON/physiq1/klev,play,JM,coef1,coef2
830      COMMON/physiq2/nblvlm,playm,hplaym
831
832      REAL psol
833      REAL val
834      INTEGER k, mlz, mlzh
835
836
837      do k=1,klev
838         val=play(k)
839       if (val .gt. playm(1)) then
840          mlz = 0
841          JM(1) = mlz
842          coef1(1)=(playm(mlz+1)-val)
843     *             /(playm(mlz+1)-psol)
844          coef2(1)=(val-psol)
845     *             /(playm(mlz+1)-psol)
846       else
847         do mlz=1,nblvlm
848          if (     val .le. playm(mlz)
849     *       .and. val .gt. playm(mlz+1))then
850           JM(k)=mlz
851           coef1(k)=(playm(mlz+1)-val)
852     *              /(playm(mlz+1)-playm(mlz))
853           coef2(k)=(val-playm(mlz))
854     *              /(playm(mlz+1)-playm(mlz))
855          endif
856c
857         enddo
858       endif
859      enddo
860c
861      if (play(klev) .le. playm(nblvlm)) then
862         mlz=nblvlm-1
863         JM(klev)=mlz
864         coef1(klev)=(playm(mlz+1)-val)
865     *            /(playm(mlz+1)-playm(mlz))
866         coef2(klev)=(val-playm(mlz))
867     *            /(playm(mlz+1)-playm(mlz))
868      endif
869c
870      print*,' '
871      print*,'         INTERPOLATION  : '
872      print*,' '
873      print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
874      print*,(JM(k),k=1,klev)
875      print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
876      print*,(JM(k),k=1,klev)
877      print*,' '
878      print*,'valeurs du premier coef d"interpolation pour les 9 niveaux
879     *: '
880      print*,(coef1(k),k=1,klev)
881      print*,' '
882      print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveau
883     *x: '
884      print*,(coef2(k),k=1,klev)
885c
886      return
887      end
888      SUBROUTINE phyredem (fichnom,dtime,radpas,co2_ppm,solaire,
889     .           rlat,rlon,tsol,tsoil,deltat,qsol,snow,
890     .           radsol,rugmer,agesno,
891     .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,
892     .           t_ancien, q_ancien)
893      IMPLICIT none
894c======================================================================
895c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
896c Objet: Ecriture de l'etat de redemarrage pour la physique
897c======================================================================
898#include "dimensions.h"
899#include "dimphy.h"
900#include "netcdf.inc"
901#include "indicesol.h"
902#include "dimsoil.h"
903#include "clesphys.h"
904#include "control.h"
905#include "temps.h"
906c======================================================================
907      CHARACTER*(*) fichnom
908      REAL dtime
909      INTEGER radpas
910      REAL rlat(klon), rlon(klon)
911      REAL co2_ppm
912      REAL solaire
913      REAL tsol(klon,nbsrf)
914      REAL tsoil(klon,nsoilmx,nbsrf)
915      REAL deltat(klon)
916      REAL qsol(klon,nbsrf)
917      REAL snow(klon,nbsrf)
918      REAL radsol(klon)
919      REAL rugmer(klon)
920      REAL agesno(klon)
921      REAL zmea(klon)
922      REAL zstd(klon)
923      REAL zsig(klon)
924      REAL zgam(klon)
925      REAL zthe(klon)
926      REAL zpic(klon)
927      REAL zval(klon)
928      REAL rugsrel(klon)
929      REAL t_ancien(klon,klev), q_ancien(klon,klev)
930c
931      INTEGER nid, nvarid, idim1, idim2, idim3
932      INTEGER ierr
933      INTEGER length
934      PARAMETER (length=100)
935      REAL tab_cntrl(length)
936c
937      INTEGER isoil, nsrf
938      CHARACTER*7 str7
939      CHARACTER*2 str2
940c
941      ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
942      IF (ierr.NE.NF_NOERR) THEN
943        write(6,*)' Pb d''ouverture du fichier '//fichnom
944        write(6,*)' ierr = ', ierr
945        CALL ABORT
946      ENDIF
947c
948      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 28,
949     .                       "Fichier redemmarage physique")
950c
951      ierr = NF_DEF_DIM (nid, "index", length, idim1)
952      ierr = NF_DEF_DIM (nid, "points_physiques", klon, idim2)
953      ierr = NF_DEF_DIM (nid, "horizon_vertical", klon*klev, idim3)
954c
955      ierr = NF_ENDDEF(nid)
956c
957      DO ierr = 1, length
958         tab_cntrl(ierr) = 0.0
959      ENDDO
960      tab_cntrl(1) = dtime
961      tab_cntrl(2) = radpas
962      tab_cntrl(3) = co2_ppm
963      tab_cntrl(4) = solaire
964      tab_cntrl(5) = iflag_con
965      tab_cntrl(6) = nbapp_rad
966
967      IF( cycle_diurne ) tab_cntrl( 7 ) = 1.
968      IF(   soil_model ) tab_cntrl( 8 ) = 1.
969      IF(     new_oliq ) tab_cntrl( 9 ) = 1.
970      IF(     ok_orodr ) tab_cntrl(10 ) = 1.
971      IF(     ok_orolf ) tab_cntrl(11 ) = 1.
972
973      tab_cntrl(13) = dayref
974      tab_cntrl(14) = anneeref
975      tab_cntrl(13) = day_end
976      tab_cntrl(14) = anne_ini
977c
978      ierr = NF_REDEF (nid)
979#ifdef NC_DOUBLE
980      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
981#else
982      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
983#endif
984      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22,
985     .                        "Parametres de controle")
986      ierr = NF_ENDDEF(nid)
987#ifdef NC_DOUBLE
988      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
989#else
990      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
991#endif
992c
993      ierr = NF_REDEF (nid)
994#ifdef NC_DOUBLE
995      ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
996#else
997      ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
998#endif
999      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 32,
1000     .                        "Longitudes de la grille physique")
1001      ierr = NF_ENDDEF(nid)
1002#ifdef NC_DOUBLE
1003      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlon)
1004#else
1005      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlon)
1006#endif
1007c
1008      ierr = NF_REDEF (nid)
1009#ifdef NC_DOUBLE
1010      ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
1011#else
1012      ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
1013#endif
1014      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 31,
1015     .                        "Latitudes de la grille physique")
1016      ierr = NF_ENDDEF(nid)
1017#ifdef NC_DOUBLE
1018      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlat)
1019#else
1020      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlat)
1021#endif
1022c
1023c
1024      DO nsrf = 1, nbsrf
1025        IF (nsrf.LE.99) THEN
1026        WRITE(str2,'(i2.2)') nsrf
1027        ierr = NF_REDEF (nid)
1028#ifdef NC_DOUBLE
1029        ierr = NF_DEF_VAR (nid, "TS"//str2, NF_DOUBLE, 1, idim2,nvarid)
1030#else
1031        ierr = NF_DEF_VAR (nid, "TS"//str2, NF_FLOAT, 1, idim2,nvarid)
1032#endif
1033        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
1034     .                        "Temperature de surface No."//str2)
1035        ierr = NF_ENDDEF(nid)
1036        ELSE
1037        PRINT*, "Trop de sous-mailles"
1038        CALL abort
1039        ENDIF
1040#ifdef NC_DOUBLE
1041        ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsol(1,nsrf))
1042#else
1043        ierr = NF_PUT_VAR_REAL (nid,nvarid,tsol(1,nsrf))
1044#endif
1045      ENDDO
1046c
1047      DO nsrf = 1, nbsrf
1048      DO isoil=1, nsoilmx
1049        IF (isoil.LE.99 .AND. nsrf.LE.99) THEN
1050        WRITE(str7,'(i2.2,"srf",i2.2)') isoil,nsrf
1051        ierr = NF_REDEF (nid)
1052#ifdef NC_DOUBLE
1053        ierr = NF_DEF_VAR (nid, "Tsoil"//str7,NF_DOUBLE,1,idim2,nvarid)
1054#else
1055        ierr = NF_DEF_VAR (nid, "Tsoil"//str7,NF_FLOAT,1,idim2,nvarid)
1056#endif
1057        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 29,
1058     .                        "Temperature du sol No."//str7)
1059        ierr = NF_ENDDEF(nid)
1060        ELSE
1061        PRINT*, "Trop de couches"
1062        CALL abort
1063        ENDIF
1064#ifdef NC_DOUBLE
1065        ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsoil(1,isoil,nsrf))
1066#else
1067        ierr = NF_PUT_VAR_REAL (nid,nvarid,tsoil(1,isoil,nsrf))
1068#endif
1069      ENDDO
1070      ENDDO
1071c
1072c
1073      ierr = NF_REDEF (nid)
1074#ifdef NC_DOUBLE
1075      ierr = NF_DEF_VAR (nid, "DELTAT", NF_DOUBLE, 1, idim2,nvarid)
1076#else
1077      ierr = NF_DEF_VAR (nid, "DELTAT", NF_FLOAT, 1, idim2,nvarid)
1078#endif
1079      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 33,
1080     .                        "Ecart de la SST (pour slab-ocean)")
1081      ierr = NF_ENDDEF(nid)
1082#ifdef NC_DOUBLE
1083      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,deltat)
1084#else
1085      ierr = NF_PUT_VAR_REAL (nid,nvarid,deltat)
1086#endif
1087c
1088      DO nsrf = 1, nbsrf
1089        IF (nsrf.LE.99) THEN
1090        WRITE(str2,'(i2.2)') nsrf
1091        ierr = NF_REDEF (nid)
1092#ifdef NC_DOUBLE
1093        ierr = NF_DEF_VAR (nid,"QS"//str2,NF_DOUBLE,1,idim2,nvarid)
1094#else
1095        ierr = NF_DEF_VAR (nid,"QS"//str2,NF_FLOAT,1,idim2,nvarid)
1096#endif
1097        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
1098     .                        "Humidite de surface No."//str2)
1099        ierr = NF_ENDDEF(nid)
1100        ELSE
1101        PRINT*, "Trop de sous-mailles"
1102        CALL abort
1103        ENDIF
1104#ifdef NC_DOUBLE
1105      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,qsol(1,nsrf))
1106#else
1107      ierr = NF_PUT_VAR_REAL (nid,nvarid,qsol(1,nsrf))
1108#endif
1109      ENDDO
1110c
1111      DO nsrf = 1, nbsrf
1112        IF (nsrf.LE.99) THEN
1113        WRITE(str2,'(i2.2)') nsrf
1114        ierr = NF_REDEF (nid)
1115#ifdef NC_DOUBLE
1116        ierr = NF_DEF_VAR (nid,"SNOW"//str2,NF_DOUBLE,1,idim2,nvarid)
1117#else
1118        ierr = NF_DEF_VAR (nid,"SNOW"//str2,NF_FLOAT,1,idim2,nvarid)
1119#endif
1120        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22,
1121     .                        "Neige de surface No."//str2)
1122        ierr = NF_ENDDEF(nid)
1123        ELSE
1124        PRINT*, "Trop de sous-mailles"
1125        CALL abort
1126        ENDIF
1127#ifdef NC_DOUBLE
1128      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,snow(1,nsrf))
1129#else
1130      ierr = NF_PUT_VAR_REAL (nid,nvarid,snow(1,nsrf))
1131#endif
1132      ENDDO
1133c
1134      ierr = NF_REDEF (nid)
1135#ifdef NC_DOUBLE
1136      ierr = NF_DEF_VAR (nid, "RADS", NF_DOUBLE, 1, idim2,nvarid)
1137#else
1138      ierr = NF_DEF_VAR (nid, "RADS", NF_FLOAT, 1, idim2,nvarid)
1139#endif
1140      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
1141     .                        "Rayonnement net a la surface")
1142      ierr = NF_ENDDEF(nid)
1143#ifdef NC_DOUBLE
1144      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,radsol)
1145#else
1146      ierr = NF_PUT_VAR_REAL (nid,nvarid,radsol)
1147#endif
1148c
1149      ierr = NF_REDEF (nid)
1150#ifdef NC_DOUBLE
1151      ierr = NF_DEF_VAR (nid, "RUGMER", NF_DOUBLE, 1, idim2,nvarid)
1152#else
1153      ierr = NF_DEF_VAR (nid, "RUGMER", NF_FLOAT, 1, idim2,nvarid)
1154#endif
1155      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
1156     .                        "Longueur de rugosite sur mer")
1157      ierr = NF_ENDDEF(nid)
1158#ifdef NC_DOUBLE
1159      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rugmer)
1160#else
1161      ierr = NF_PUT_VAR_REAL (nid,nvarid,rugmer)
1162#endif
1163c
1164      ierr = NF_REDEF (nid)
1165#ifdef NC_DOUBLE
1166      ierr = NF_DEF_VAR (nid, "AGESNO", NF_DOUBLE, 1, idim2,nvarid)
1167#else
1168      ierr = NF_DEF_VAR (nid, "AGESNO", NF_FLOAT, 1, idim2,nvarid)
1169#endif
1170      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 15,
1171     .                        "Age de la neige")
1172      ierr = NF_ENDDEF(nid)
1173#ifdef NC_DOUBLE
1174      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,agesno)
1175#else
1176      ierr = NF_PUT_VAR_REAL (nid,nvarid,agesno)
1177#endif
1178c
1179      ierr = NF_REDEF (nid)
1180#ifdef NC_DOUBLE
1181      ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
1182#else
1183      ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
1184#endif
1185      ierr = NF_ENDDEF(nid)
1186#ifdef NC_DOUBLE
1187      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmea)
1188#else
1189      ierr = NF_PUT_VAR_REAL (nid,nvarid,zmea)
1190#endif
1191c
1192      ierr = NF_REDEF (nid)
1193#ifdef NC_DOUBLE
1194      ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
1195#else
1196      ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
1197#endif
1198      ierr = NF_ENDDEF(nid)
1199#ifdef NC_DOUBLE
1200      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zstd)
1201#else
1202      ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
1203#endif
1204c
1205      ierr = NF_REDEF (nid)
1206#ifdef NC_DOUBLE
1207      ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
1208#else
1209      ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
1210#endif
1211      ierr = NF_ENDDEF(nid)
1212#ifdef NC_DOUBLE
1213      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zsig)
1214#else
1215      ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
1216#endif
1217c
1218      ierr = NF_REDEF (nid)
1219#ifdef NC_DOUBLE
1220      ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
1221#else
1222      ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
1223#endif
1224      ierr = NF_ENDDEF(nid)
1225#ifdef NC_DOUBLE
1226      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zgam)
1227#else
1228      ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
1229#endif
1230c
1231      ierr = NF_REDEF (nid)
1232#ifdef NC_DOUBLE
1233      ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
1234#else
1235      ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
1236#endif
1237      ierr = NF_ENDDEF(nid)
1238#ifdef NC_DOUBLE
1239      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zthe)
1240#else
1241      ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
1242#endif
1243c
1244      ierr = NF_REDEF (nid)
1245#ifdef NC_DOUBLE
1246      ierr = NF_DEF_VAR (nid, "ZPIC", NF_DOUBLE, 1, idim2,nvarid)
1247#else
1248      ierr = NF_DEF_VAR (nid, "ZPIC", NF_FLOAT, 1, idim2,nvarid)
1249#endif
1250      ierr = NF_ENDDEF(nid)
1251#ifdef NC_DOUBLE
1252      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zpic)
1253#else
1254      ierr = NF_PUT_VAR_REAL (nid,nvarid,zpic)
1255#endif
1256c
1257      ierr = NF_REDEF (nid)
1258#ifdef NC_DOUBLE
1259      ierr = NF_DEF_VAR (nid, "ZVAL", NF_DOUBLE, 1, idim2,nvarid)
1260#else
1261      ierr = NF_DEF_VAR (nid, "ZVAL", NF_FLOAT, 1, idim2,nvarid)
1262#endif
1263      ierr = NF_ENDDEF(nid)
1264#ifdef NC_DOUBLE
1265      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zval)
1266#else
1267      ierr = NF_PUT_VAR_REAL (nid,nvarid,zval)
1268#endif
1269c
1270      ierr = NF_REDEF (nid)
1271#ifdef NC_DOUBLE
1272      ierr = NF_DEF_VAR (nid, "RUGSREL", NF_DOUBLE, 1, idim2,nvarid)
1273#else
1274      ierr = NF_DEF_VAR (nid, "RUGSREL", NF_FLOAT, 1, idim2,nvarid)
1275#endif
1276      ierr = NF_ENDDEF(nid)
1277#ifdef NC_DOUBLE
1278      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rugsrel)
1279#else
1280      ierr = NF_PUT_VAR_REAL (nid,nvarid,rugsrel)
1281#endif
1282c
1283      ierr = NF_REDEF (nid)
1284#ifdef NC_DOUBLE
1285      ierr = NF_DEF_VAR (nid, "TANCIEN", NF_DOUBLE, 1, idim3,nvarid)
1286#else
1287      ierr = NF_DEF_VAR (nid, "TANCIEN", NF_FLOAT, 1, idim3,nvarid)
1288#endif
1289      ierr = NF_ENDDEF(nid)
1290#ifdef NC_DOUBLE
1291      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,t_ancien)
1292#else
1293      ierr = NF_PUT_VAR_REAL (nid,nvarid,t_ancien)
1294#endif
1295c
1296      ierr = NF_REDEF (nid)
1297#ifdef NC_DOUBLE
1298      ierr = NF_DEF_VAR (nid, "QANCIEN", NF_DOUBLE, 1, idim3,nvarid)
1299#else
1300      ierr = NF_DEF_VAR (nid, "QANCIEN", NF_FLOAT, 1, idim3,nvarid)
1301#endif
1302      ierr = NF_ENDDEF(nid)
1303#ifdef NC_DOUBLE
1304      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q_ancien)
1305#else
1306      ierr = NF_PUT_VAR_REAL (nid,nvarid,q_ancien)
1307#endif
1308c
1309      ierr = NF_CLOSE(nid)
1310c
1311      RETURN
1312      END
1313      subroutine physdem(lonfi, latfi,phystep,radpas,co2_ppm,
1314     .                   solaire, ts, ws, 
1315     .                   sn, radsol, deltat, rugmer,
1316     .                   agesno, zmea, zstd, zsig,
1317     .                   zgam, zthe, zpic, zval,
1318     .                   rugsrel)
1319
1320      IMPLICIT none
1321c-------------------------------------------------------------
1322C Author : L. Fairhead
1323C Date   : 01/10/1999
1324C Objet  : Ecriture des etats initiaux physiques
1325c-------------------------------------------------------------
1326c
1327c
1328c
1329      INTEGER ivap
1330      PARAMETER (ivap=1)
1331c
1332      REAL qsolmax
1333      PARAMETER ( qsolmax = 150.0 )
1334c
1335#include "dimensions.h"
1336#include "paramet.h"
1337#include "dimphy.h"
1338#include "control.h"
1339#include "netcdf.inc"
1340c
1341      INTEGER nid
1342
1343c Ajout de quelques parametres orographiques (F. LOTT janvier 1995)
1344
1345      REAL zmea(iip1,jjp1),zstd(iip1,jjp1)
1346      REAL zsig(iip1,jjp1),zgam(iip1,jjp1),zthe(iip1,jjp1)
1347      REAL zpic(iip1,jjp1),zval(iip1,jjp1)
1348      REAL rugsrel(iip1,jjp1)
1349      INTEGER idayref,anneeref
1350
1351
1352      integer ierr, idim1, idim2, nvarid
1353
1354c
1355      REAL phystep
1356      INTEGER radpas
1357      REAL co2_ppm
1358      REAL solaire
1359      REAL latfi(klon), lonfi(klon)
1360      REAL champhys(klon)
1361      REAL ts(klon)
1362      REAL deltat(klon)
1363      REAL ws(klon)
1364      REAL sn(klon)
1365      REAL radsol(klon)
1366      REAL rugmer(klon)
1367      REAL agesno(klon)
1368      INTEGER length
1369      PARAMETER (length=100)
1370      REAL tab_cntrl(length)
1371      real pi
1372
1373c
1374
1375#include "serre.h"
1376#include "clesphys.h"
1377#include "fxyprim.h"
1378c-----------------------------------------------------------------------
1379c
1380c  stockage sur le fichier Physique:
1381c
1382      pi=2.*asin(1.)
1383      ierr = NF_CREATE("startphy.nc", NF_CLOBBER, nid)
1384      IF (ierr.NE.NF_NOERR) THEN
1385        WRITE(6,*)' Pb d''ouverture du fichier startphy.nc'
1386        WRITE(6,*)' ierr = ', ierr
1387        CALL ABORT
1388      ENDIF
1389c
1390      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 28,
1391     .                       "Fichier demmarage physique")
1392c
1393      ierr = NF_DEF_DIM (nid, "index", length, idim1)
1394      ierr = NF_DEF_DIM (nid, "points_physiques", klon, idim2)
1395c
1396      ierr = NF_ENDDEF(nid)
1397c
1398      DO ierr = 1, length
1399         tab_cntrl(ierr) = 0.0
1400      ENDDO
1401      tab_cntrl(1)  = phystep
1402      tab_cntrl(2)  = radpas
1403      tab_cntrl(3)  = co2_ppm
1404      tab_cntrl(4)  = solaire
1405      tab_cntrl(5)  = iflag_con
1406      tab_cntrl(6)  = nbapp_rad
1407c
1408cc     Modif ( P. Le Van )
1409c
1410       tab_cntrl( 7 ) = 0. 
1411       tab_cntrl( 8 ) = 0. 
1412       tab_cntrl( 9 ) = 0. 
1413       tab_cntrl(10 ) = 0. 
1414       tab_cntrl(11 ) = 0. 
1415       tab_cntrl(12 ) = 0. 
1416
1417      IF(  cycle_diurne )  tab_cntrl( 7 ) = 1. 
1418      IF(   soil_model  )  tab_cntrl( 8 ) = 1. 
1419      IF(    new_oliq   )  tab_cntrl( 9 ) = 1. 
1420      IF(    ok_orodr   )  tab_cntrl(10 ) = 1. 
1421      IF(    ok_orolf   )  tab_cntrl(11 ) = 1. 
1422      IF(  ok_limitvrai )  tab_cntrl(12 ) = 1. 
1423
1424      tab_cntrl(13)  = dayref
1425      tab_cntrl(14)  = anneeref
1426
1427
1428cc   ***    new_oliq   (  commentaires de L. LI dans routine physique ) 
1429cc   ***  ok_orodr  et ok_orolf   si on appelle l'orographie      ****
1430
1431c
1432      ierr = NF_REDEF (nid)
1433#ifdef NC_DOUBLE
1434      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
1435#else
1436      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
1437#endif
1438      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22,
1439     .                        "Parametres de controle")
1440      ierr = NF_ENDDEF(nid)
1441#ifdef NC_DOUBLE
1442      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
1443#else
1444      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
1445#endif
1446c
1447      ierr = NF_REDEF (nid)
1448#ifdef NC_DOUBLE
1449      ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
1450#else
1451      ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
1452#endif
1453      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 32,
1454     .                        "Longitudes de la grille physique")
1455      ierr = NF_ENDDEF(nid)
1456
1457#ifdef NC_DOUBLE
1458      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lonfi)
1459#else
1460      ierr = NF_PUT_VAR_REAL (nid,nvarid,lonfi)
1461#endif
1462c
1463      ierr = NF_REDEF (nid)
1464#ifdef NC_DOUBLE
1465      ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
1466#else
1467      ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
1468#endif
1469      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 31,
1470     .                        "Latitudes de la grille physique")
1471      ierr = NF_ENDDEF(nid)
1472#ifdef NC_DOUBLE
1473      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,latfi)
1474#else
1475      ierr = NF_PUT_VAR_REAL (nid,nvarid,latfi)
1476#endif
1477c
1478      ierr = NF_REDEF (nid)
1479#ifdef NC_DOUBLE
1480      ierr = NF_DEF_VAR (nid, "TS", NF_DOUBLE, 1, idim2,nvarid)
1481#else
1482      ierr = NF_DEF_VAR (nid, "TS", NF_FLOAT, 1, idim2,nvarid)
1483#endif
1484      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
1485     .                        "Temperature de la surface")
1486      ierr = NF_ENDDEF(nid)
1487#ifdef NC_DOUBLE
1488      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ts)
1489#else
1490      ierr = NF_PUT_VAR_REAL (nid,nvarid,ts)
1491#endif
1492c
1493      ierr = NF_REDEF (nid)
1494#ifdef NC_DOUBLE
1495      ierr = NF_DEF_VAR (nid, "QS", NF_DOUBLE, 1, idim2,nvarid)
1496#else
1497      ierr = NF_DEF_VAR (nid, "QS", NF_FLOAT, 1, idim2,nvarid)
1498#endif
1499      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 15,
1500     .                        "Humidite du sol")
1501      ierr = NF_ENDDEF(nid)
1502#ifdef NC_DOUBLE
1503      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ws)
1504#else
1505      ierr = NF_PUT_VAR_REAL (nid,nvarid,ws)
1506#endif
1507c
1508      ierr = NF_REDEF (nid)
1509#ifdef NC_DOUBLE
1510      ierr = NF_DEF_VAR (nid, "SNOW", NF_DOUBLE, 1, idim2,nvarid)
1511#else
1512      ierr = NF_DEF_VAR (nid, "SNOW", NF_FLOAT, 1, idim2,nvarid)
1513#endif
1514      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 5,
1515     .                        "Neige")
1516      ierr = NF_ENDDEF(nid)
1517#ifdef NC_DOUBLE
1518      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,sn)
1519#else
1520      ierr = NF_PUT_VAR_REAL (nid,nvarid,sn)
1521#endif
1522c
1523      ierr = NF_REDEF (nid)
1524#ifdef NC_DOUBLE
1525      ierr = NF_DEF_VAR (nid, "RADS", NF_DOUBLE, 1, idim2,nvarid)
1526#else
1527      ierr = NF_DEF_VAR (nid, "RADS", NF_FLOAT, 1, idim2,nvarid)
1528#endif
1529      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
1530     .                        "Rayonnement net a la surface")
1531      ierr = NF_ENDDEF(nid)
1532#ifdef NC_DOUBLE
1533      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,radsol)
1534#else
1535      ierr = NF_PUT_VAR_REAL (nid,nvarid,radsol)
1536#endif
1537c
1538      ierr = NF_REDEF (nid)
1539#ifdef NC_DOUBLE
1540      ierr = NF_DEF_VAR (nid, "DELTAT", NF_DOUBLE, 1, idim2,nvarid)
1541#else
1542      ierr = NF_DEF_VAR (nid, "DELTAT", NF_FLOAT, 1, idim2,nvarid)
1543#endif
1544      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 33,
1545     .                        "Ecart de la SST (pour slab-ocean)")
1546      ierr = NF_ENDDEF(nid)
1547#ifdef NC_DOUBLE
1548      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,deltat)
1549#else
1550      ierr = NF_PUT_VAR_REAL (nid,nvarid,deltat)
1551#endif
1552c
1553      ierr = NF_REDEF (nid)
1554#ifdef NC_DOUBLE
1555      ierr = NF_DEF_VAR (nid, "RUGMER", NF_DOUBLE, 1, idim2,nvarid)
1556#else
1557      ierr = NF_DEF_VAR (nid, "RUGMER", NF_FLOAT, 1, idim2,nvarid)
1558#endif
1559      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
1560     .                        "Longueur de rugosite sur mer")
1561      ierr = NF_ENDDEF(nid)
1562#ifdef NC_DOUBLE
1563      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rugmer)
1564#else
1565      ierr = NF_PUT_VAR_REAL (nid,nvarid,rugmer)
1566#endif
1567c
1568      ierr = NF_REDEF (nid)
1569#ifdef NC_DOUBLE
1570      ierr = NF_DEF_VAR (nid, "AGESNO", NF_DOUBLE, 1, idim2,nvarid)
1571#else
1572      ierr = NF_DEF_VAR (nid, "AGESNO", NF_FLOAT, 1, idim2,nvarid)
1573#endif
1574      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 15,
1575     .                        "Age de la neige")
1576      ierr = NF_ENDDEF(nid)
1577#ifdef NC_DOUBLE
1578      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,agesno)
1579#else
1580      ierr = NF_PUT_VAR_REAL (nid,nvarid,agesno)
1581#endif
1582c
1583      CALL gr_dyn_fi(1, iip1, jjp1, klon, zmea, champhys)
1584      ierr = NF_REDEF (nid)
1585#ifdef NC_DOUBLE
1586      ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
1587#else
1588      ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
1589#endif
1590      ierr = NF_ENDDEF(nid)
1591#ifdef NC_DOUBLE
1592      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
1593#else
1594      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
1595#endif
1596c
1597      CALL gr_dyn_fi(1, iip1, jjp1, klon, zstd, champhys)
1598      ierr = NF_REDEF (nid)
1599#ifdef NC_DOUBLE
1600      ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
1601#else
1602      ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
1603#endif
1604      ierr = NF_ENDDEF(nid)
1605#ifdef NC_DOUBLE
1606      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
1607#else
1608      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
1609#endif
1610
1611      CALL gr_dyn_fi(1, iip1, jjp1, klon, zsig, champhys)
1612      ierr = NF_REDEF (nid)
1613#ifdef NC_DOUBLE
1614      ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
1615#else
1616      ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
1617#endif
1618      ierr = NF_ENDDEF(nid)
1619#ifdef NC_DOUBLE
1620      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
1621#else
1622      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
1623#endif
1624
1625      CALL gr_dyn_fi(1, iip1, jjp1, klon, zgam, champhys)
1626      ierr = NF_REDEF (nid)
1627#ifdef NC_DOUBLE
1628      ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
1629#else
1630      ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
1631#endif
1632      ierr = NF_ENDDEF(nid)
1633#ifdef NC_DOUBLE
1634      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
1635#else
1636      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
1637#endif
1638
1639      CALL gr_dyn_fi(1, iip1, jjp1, klon, zthe, champhys)
1640      ierr = NF_REDEF (nid)
1641#ifdef NC_DOUBLE
1642      ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
1643#else
1644      ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
1645#endif
1646      ierr = NF_ENDDEF(nid)
1647#ifdef NC_DOUBLE
1648      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
1649#else
1650      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
1651#endif
1652
1653      CALL gr_dyn_fi(1, iip1, jjp1, klon, zpic, champhys)
1654      ierr = NF_REDEF (nid)
1655#ifdef NC_DOUBLE
1656      ierr = NF_DEF_VAR (nid, "ZPIC", NF_DOUBLE, 1, idim2,nvarid)
1657#else
1658      ierr = NF_DEF_VAR (nid, "ZPIC", NF_FLOAT, 1, idim2,nvarid)
1659#endif
1660      ierr = NF_ENDDEF(nid)
1661#ifdef NC_DOUBLE
1662      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
1663#else
1664      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
1665#endif
1666
1667      CALL gr_dyn_fi(1, iip1, jjp1, klon, zval, champhys)
1668      ierr = NF_REDEF (nid)
1669#ifdef NC_DOUBLE
1670      ierr = NF_DEF_VAR (nid, "ZVAL", NF_DOUBLE, 1, idim2,nvarid)
1671#else
1672      ierr = NF_DEF_VAR (nid, "ZVAL", NF_FLOAT, 1, idim2,nvarid)
1673#endif
1674      ierr = NF_ENDDEF(nid)
1675#ifdef NC_DOUBLE
1676      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
1677#else
1678      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
1679#endif
1680
1681      CALL gr_dyn_fi(1, iip1, jjp1, klon, rugsrel, champhys)
1682      ierr = NF_REDEF (nid)
1683#ifdef NC_DOUBLE
1684      ierr = NF_DEF_VAR (nid, "RUGSREL", NF_DOUBLE, 1, idim2,nvarid)
1685#else
1686      ierr = NF_DEF_VAR (nid, "RUGSREL", NF_FLOAT, 1, idim2,nvarid)
1687#endif
1688      ierr = NF_ENDDEF(nid)
1689#ifdef NC_DOUBLE
1690      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
1691#else
1692      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
1693#endif
1694c
1695      ierr = NF_CLOSE(nid)
1696
1697      RETURN
1698
1699      END
1700*CMZ :          28/02/95  17.58.56  by  Unknown
1701*-- Author :
1702      CHARACTER*(*) FUNCTION SPACES(STR,NSPACE)
1703C
1704C CERN PROGLIB# M433    SPACES          .VERSION KERNFOR  4.14  860211
1705C ORIG.  6/05/86 M.GOOSSENS/DD
1706C
1707C-    The function value SPACES returns the character string STR with
1708C-    leading blanks removed and each occurence of one or more blanks
1709C-    replaced by NSPACE blanks inside the string STR
1710C
1711      CHARACTER*(*) STR
1712C
1713      LENSPA = LEN(SPACES)
1714      SPACES = ' '
1715      IF (NSPACE.LT.0) NSPACE = 0
1716      IBLANK = 1
1717      ISPACE = 1
1718  100 INONBL = INDEXC(STR(IBLANK:),' ')
1719      IF (INONBL.EQ.0) THEN
1720          SPACES(ISPACE:) = STR(IBLANK:)
1721                                                    GO TO 999
1722      ENDIF
1723      INONBL = INONBL + IBLANK - 1
1724      IBLANK = INDEX(STR(INONBL:),' ')
1725      IF (IBLANK.EQ.0) THEN
1726          SPACES(ISPACE:) = STR(INONBL:)
1727                                                    GO TO 999
1728      ENDIF
1729      IBLANK = IBLANK + INONBL - 1
1730      SPACES(ISPACE:) = STR(INONBL:IBLANK-1)
1731      ISPACE = ISPACE + IBLANK - INONBL + NSPACE
1732      IF (ISPACE.LE.LENSPA)                         GO TO 100
1733  999 END
1734      FUNCTION INDEXC(STR,SSTR)
1735C
1736C CERN PROGLIB# M433    INDEXC          .VERSION KERNFOR  4.14  860211
1737C ORIG. 26/03/86 M.GOOSSENS/DD
1738C
1739C-    Find the leftmost position where substring SSTR does not match
1740C-    string STR scanning forward
1741C
1742      CHARACTER*(*) STR,SSTR
1743C
1744      LENS   = LEN(STR)
1745      LENSS  = LEN(SSTR)
1746C
1747      DO 10 I=1,LENS-LENSS+1
1748          IF (STR(I:I+LENSS-1).NE.SSTR) THEN
1749              INDEXC = I
1750                                         GO TO 999
1751          ENDIF
1752   10 CONTINUE
1753      INDEXC = 0
1754C
1755  999 END
Note: See TracBrowser for help on using the repository browser.