source: LMDZ.3.3/trunk/libf/dyn3d/nudge.F @ 2231

Last change on this file since 2231 was 356, checked in by lmdz, 23 years ago

Lecture du fichier netcdf pour determination du nombre de niveaux verticaux MAF
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 7.6 KB
Line 
1c $Header$
2      subroutine nudge(itau,ucov,vcov,teta,masse,ps)
3
4      IMPLICIT NONE
5
6c      ......   Version  du 10/01/98    ..........
7
8c             avec  coordonnees  verticales hybrides
9c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
10
11c=======================================================================
12c
13c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
14c   -------
15c
16c   Objet:
17c   ------
18c
19c   GCM LMD nouvelle grille
20c
21c=======================================================================
22
23c   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv
24c        et possibilite d'appeler une fonction f(y)  a derivee tangente
25c        hyperbolique a la  place de la fonction a derivee sinusoidale.         
26
27c   ...  Possibilite de choisir le shema de Van-leer pour l'advection de
28c         q  , en faisant iadv = 10  dans   traceur  (29/04/97) .
29c
30c-----------------------------------------------------------------------
31c   Declarations:
32c   -------------
33
34#include "dimensions.h"
35#include "paramet.h"
36#include "comconst.h"
37#include "comdissnew.h"
38#include "comvert.h"
39#include "comgeom.h"
40#include "logic.h"
41#include "temps.h"
42#include "control.h"
43#include "ener.h"
44#include "netcdf.inc"
45#include "description.h"
46#include "serre.h"
47#include "tracstoke.h"
48
49
50c   variables dynamiques
51      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
52      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
53      REAL ps(ip1jmp1)                       ! pression  au sol
54      REAL masse(ip1jmp1,llm)                ! masse d'air
55
56c   variables dynamiques pour les reanalyses.
57      REAL ucovrea1(ip1jmp1,llm),vcovrea1(ip1jm,llm) !vts cov reas
58      REAL tetarea1(ip1jmp1,llm)             ! temp pot  reales
59      REAL masserea1(ip1jmp1,llm)             ! masse
60      REAL psrea1(ip1jmp1)             ! ps
61      REAL ucovrea2(ip1jmp1,llm),vcovrea2(ip1jm,llm) !vts cov reas
62      REAL tetarea2(ip1jmp1,llm)             ! temp pot  reales
63      REAL masserea2(ip1jmp1,llm)             ! masse
64      REAL psrea2(ip1jmp1)             ! ps
65      real alphaT,alphaP,alphau,alphav
66      real dday_step,toto,reste,itau_test
67      INTEGER step_rea,count_no_rea
68
69
70c
71      INTEGER itau,ij,l
72      integer ncidt,varidpl,nlev,status
73      integer rcod,rid
74      real ditau,tau,a
75      save nlev
76
77      LOGICAL first,ncep
78      integer online
79      save first,online,ncep
80      data first,online/.true.,1/
81      data ncep/.false./
82
83      save ucovrea1,vcovrea1,tetarea1,masserea1,psrea1
84      save ucovrea2,vcovrea2,tetarea2,masserea2,psrea2
85
86      save alphaT,alphau,alphav,alphaP,itau_test
87      save step_rea,count_no_rea
88
89c-----------------------------------------------------------------------
90c   initialisations pour la lecture des reanalyses.
91c    alpha determine la part des injections de donnees a chaque etape
92c    alpha=1 signifie pas d'injection
93c    alpha=0 signifie injection totale
94c-----------------------------------------------------------------------
95
96      print*,'ONLINE=',online
97      if(online.eq.-1) then
98          return
99      endif
100
101      if (first) then
102         print*
103     s   ,'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
104cnec         read(*,*) online
105cnec      online=-1
106cnec     print*,'Entrer les constantes de temps de rappel en jours'
107cnec         read(*,*)alphaT
108cnec         read(*,*)alphau
109cnec         read(*,*)alphav
110cnec         read(*,*)alphaP
111             alphaT=0.1
112             alphau=0.1
113             alphav=0.1
114             alphaP=1.e10
115         print*,'alphaT,alphau,alphav,alphaP'
116     s          ,alphaT,alphau,alphav,alphaP
117         if(online.eq.-1) return
118         print*,'alpha rappel pour T, u, v, P ',
119     s    alphaT,alphau,alphav,alphaP
120c   En fait, le coef alpha qu'on utilise est
121c   x_gcm = a * x_gcm + ( 1 - a ) * x_reanalys
122c   on a alors a=1-dt/tau
123c   ou dt est la frequence a laquelle on corrige (ici iperiod*dtvr)
124c   et tau la constante de temps lue. D'ou:
125         if (online.eq.1) then
126            alphaT=1.-iperiod*dtvr/(daysec*alphaT)
127            alphau=1.-iperiod*dtvr/(daysec*alphau)
128            alphav=1.-iperiod*dtvr/(daysec*alphav)
129            alphaP=1.-iperiod*dtvr/(daysec*alphaP)
130         else
131            alphaT=0.
132            alphau=0.
133            alphav=0.
134            alphaP=0.
135c           physic=.false.
136         endif
137         print*,'alpha rappel pour T, u, v, P ',
138     s    alphaT,alphau,alphav,alphaP
139
140         itau_test=1001
141         step_rea=1
142         count_no_rea=0
143
144c    itau_test    montre si l'importation a deja ete faite au rang itau
145c lecture d'un fichier netcdf pour determiner le nombre de niveaux
146         ncidt=NCOPN('T.nc',NCNOWRIT,rcod)
147         if (ncep) then
148          status=NF_INQ_DIMID(ncidt,'LEVEL',rid)
149         else
150          status=NF_INQ_DIMID(ncidt,'PRESSURE',rid)
151         endif
152          status=NF_INQ_DIMLEN(ncidt,rid,nlev)
153         print *,'nlev', nlev
154          call ncclos(ncidt,rcod)
155c   Lecture du premier etat des reanalyses.
156         call read_reanalyse(1
157     s   ,ucovrea2,vcovrea2,tetarea2,masserea2,ps,1,nlev)
158
159c-----------------------------------------------------------------------
160c   Debut de l'integration temporelle:
161c   ----------------------------------
162
163         first=.false.
164      endif ! first
165c
166C-----------------------------------------------------------------------
167C----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:
168C-----------------------------------------------------------------------
169
170      ditau=real(itau)
171      dday_step=real(day_step)
172      write(*,*)'ditau,dday_step'
173      write(*,*)ditau,dday_step
174      toto=4*ditau/dday_step
175      reste=toto-aint(toto)
176c     write(*,*)'toto,reste',toto,reste
177
178      if (reste.eq.0.) then
179        if (itau_test.eq.itau) then
180          write(*,*)'deuxieme passage de advreel a itau=',itau
181          stop
182        else
183           CALL SCOPY( ijmllm ,vcovrea2, 1, vcovrea1 , 1 )
184           CALL SCOPY( ijp1llm,ucovrea2, 1, ucovrea1 , 1 )
185           CALL SCOPY( ijp1llm,tetarea2,1,tetarea1   , 1 )
186c          CALL SCOPY( ijp1llm,masserea2,1,masserea1   , 1 )
187c          CALL SCOPY( ip1jmp1,psrea2, 1,  psrea1, 1 )
188
189          print*,'LECTURE REANALYSES, pas ',step_rea
190     s         ,'apres ',count_no_rea,' non lectures'
191           step_rea=step_rea+1
192           itau_test=itau
193           call read_reanalyse(step_rea
194     s     ,ucovrea2,vcovrea2,tetarea2,masserea2,ps,1,nlev)
195        endif
196      else
197        count_no_rea=count_no_rea+1
198      endif
199 
200C-----------------------------------------------------------------------
201c   Corrrections
202c
203c    x_gcm = a * x_gcm + (1-a) * x_reanalyses
204C-----------------------------------------------------------------------
205
206      ditau=real(itau)
207      dday_step=real(day_step)
208
209c     print*,'alpha1',alphaT,alphaP,alphau,alphav
210
211      tau=4*ditau/dday_step
212      tau=tau-aint(tau)
213
214      do l=1,llm
215         do ij=1,ip1jmp1
216            a=(1.-tau)*ucovrea1(ij,l)+tau*ucovrea2(ij,l)
217            ucov(ij,l)=alphau*ucov(ij,l)+(1-alphau)*a
218            a=(1.-tau)*tetarea1(ij,l)+tau*tetarea2(ij,l)
219            teta(ij,l)=alphaT*teta(ij,l)+(1-alphaT)*a
220c           a=(1.-tau)*masserea1(ij,l)+tau*masserea2(ij,l)
221c           masse(ij,l)=alphaP*masse(ij,l)+(1-alphaP)*a
222         enddo
223         do ij=1,ip1jm
224            a=(1.-tau)*vcovrea1(ij,l)+tau*vcovrea2(ij,l)
225            vcov(ij,l)=alphav*vcov(ij,l)+(1-alphav)*a
226         enddo
227      enddo
228
229c     call dump2d(iip1,jjp1,tetarea1,'TETA REA 1     ')
230c     call dump2d(iip1,jjp1,tetarea2,'TETA REA 2     ')
231c     call dump2d(iip1,jjp1,teta,'TETA           ')
232Cmaf  on ne nudge pas sur la pression de surface
233c      do ij=1,ip1jmp1
234c         a=(1.-tau)*psrea1(ij)+tau*psrea2(ij)
235c         ps(ij)=alphaP*ps(ij)+(1-alphaP)*a
236c      enddo
237
238
239      return
240      end
Note: See TracBrowser for help on using the repository browser.