| 1 | PROGRAM anldoppler2 |
|---|
| 2 | IMPLICIT NONE |
|---|
| 3 | c====================================================================== |
|---|
| 4 | c |
|---|
| 5 | c |
|---|
| 6 | c Programme pour analyser vent a heure locale fixee |
|---|
| 7 | c a partir d'un fichier Netcdf de la MCD |
|---|
| 8 | c |
|---|
| 9 | c======================================================================= |
|---|
| 10 | c----------------------------------------------------------------------- |
|---|
| 11 | c declarations: |
|---|
| 12 | c ------------- |
|---|
| 13 | |
|---|
| 14 | |
|---|
| 15 | #include "dimensions.h" |
|---|
| 16 | #include "paramet.h" |
|---|
| 17 | #include "comconst.h" |
|---|
| 18 | #include "comdissip.h" |
|---|
| 19 | #include "comvert.h" |
|---|
| 20 | #include "comgeom2.h" |
|---|
| 21 | #include "logic.h" |
|---|
| 22 | #include "temps.h" |
|---|
| 23 | #include "control.h" |
|---|
| 24 | #include "ener.h" |
|---|
| 25 | #include "netcdf.inc" |
|---|
| 26 | #include "description.h" |
|---|
| 27 | |
|---|
| 28 | |
|---|
| 29 | INTEGER point |
|---|
| 30 | PARAMETER(point=11) ! nbre de points d''observation |
|---|
| 31 | real olon(point),olat(point) ! Coordonnees des points observes |
|---|
| 32 | INTEGER iraie ! pour le choix de la raie du CO |
|---|
| 33 | |
|---|
| 34 | c-------------------------------------------------------- |
|---|
| 35 | c coordonees lon/lat des points observes |
|---|
| 36 | c-------------------------------------------------------- |
|---|
| 37 | |
|---|
| 38 | c ANNEE 2001 |
|---|
| 39 | |
|---|
| 40 | DATA olon/0,-30,30,0,-90,90,-75,75,0,-70,70/ |
|---|
| 41 | DATA olat/-1,0,0,40,0,0,40,40,90,-40,40/ |
|---|
| 42 | |
|---|
| 43 | c ANNEE 1999 |
|---|
| 44 | |
|---|
| 45 | c DATA olon/0,0,-30,30,-85,85,0,-90,90,0,-65,65,0/ |
|---|
| 46 | c DATA olat/18,-10,0,0,0,0,55,40,40,90,-40,-40,-70/ |
|---|
| 47 | |
|---|
| 48 | c ANNEE 1997 |
|---|
| 49 | |
|---|
| 50 | c DATA olon/0,-85,85,0,0/ |
|---|
| 51 | c DATA olat/24,0,0,90,-70/ |
|---|
| 52 | |
|---|
| 53 | c ANNEE 1992 |
|---|
| 54 | |
|---|
| 55 | c DATA olon/0,-85,85,0,0,-75,75,-90,90,0,0/ |
|---|
| 56 | c DATA olat/8,0,0,-40,40,-40,-40,40,40,90,-90/ |
|---|
| 57 | |
|---|
| 58 | c ANNEE 1990 |
|---|
| 59 | |
|---|
| 60 | c DATA olon/0,-90,90,0,0/ |
|---|
| 61 | c DATA olat/-1,0,0,70,-90/ |
|---|
| 62 | |
|---|
| 63 | c ANNEE 1988 |
|---|
| 64 | |
|---|
| 65 | c DATA olon/0,-90,90,0,0,0,-35,35,0,0,-45,45,-60,60/ |
|---|
| 66 | c DATA olat/-21,0,0,-90,60,0,0,0,-35,25,-35,-35,25,25/ |
|---|
| 67 | |
|---|
| 68 | c---------------------------------------------------------- |
|---|
| 69 | |
|---|
| 70 | INTEGER mcdreadnc |
|---|
| 71 | external mcdreadnc |
|---|
| 72 | INTEGER itau,nbpas,nbpasmx |
|---|
| 73 | PARAMETER(nbpasmx=1000000) |
|---|
| 74 | REAL temps(nbpasmx) |
|---|
| 75 | INTEGER unitlec |
|---|
| 76 | INTEGER i,j,l,itautot |
|---|
| 77 | |
|---|
| 78 | c Declarations DRS: |
|---|
| 79 | c ----------------- |
|---|
| 80 | INTEGER ierr, dimtype |
|---|
| 81 | CHARACTER*120 dimsou, dimnam*16, dimtit*80 |
|---|
| 82 | CHARACTER*40 dimunit |
|---|
| 83 | CHARACTER*100 varname |
|---|
| 84 | |
|---|
| 85 | c variables meteo dans la grille verticale GCM |
|---|
| 86 | c -------------------------------------------- |
|---|
| 87 | REAL v(iip1,jjp1,llm),u(iip1,jjp1,llm) |
|---|
| 88 | REAL t(iip1,jjp1,llm) |
|---|
| 89 | REAL ps(iip1,jjp1) , tsurf(iip1,jjp1) |
|---|
| 90 | |
|---|
| 91 | c variables meteo en coordonnee de pression |
|---|
| 92 | c -------------------------------------------- |
|---|
| 93 | REAL vp(iip1,jjp1,llm),up(iip1,jjp1,llm) |
|---|
| 94 | REAL tp(iip1,jjp1,llm) |
|---|
| 95 | |
|---|
| 96 | c Niveaux de pression (Pa) |
|---|
| 97 | real pref(llm) |
|---|
| 98 | |
|---|
| 99 | c version 25 couches |
|---|
| 100 | |
|---|
| 101 | c data pref/1291.37,1005.720,783.255,600.,475.069,350., |
|---|
| 102 | c & 288.144,174.768,106.002,64.2935,38.9960, |
|---|
| 103 | c & 23.6523,14.3458,8.70000,5.28000,3.20000, |
|---|
| 104 | c & 1.94000,1.18000,0.714000,0.433000,0.263000, |
|---|
| 105 | c & 0.118000,0.0614000,0.0333000,0.0163000/ |
|---|
| 106 | |
|---|
| 107 | c version 32 couches |
|---|
| 108 | |
|---|
| 109 | data pref/1291.37,1005.720,783.255,600.,475.069,350., |
|---|
| 110 | & 288.144,174.768,106.002,64.2935,38.9960, |
|---|
| 111 | & 23.6523,14.3458,8.70000,5.28000,3.20000, |
|---|
| 112 | & 1.94000,1.18000,0.714000,0.433000,0.263000, |
|---|
| 113 | & 0.118000,0.0614000,0.0333000,0.0163000, |
|---|
| 114 | & 8.e-3,4.e-3,2.e-3,1.e-3,5.e-4,2.e-4,1.e-4/ |
|---|
| 115 | |
|---|
| 116 | c Local: |
|---|
| 117 | c ------ |
|---|
| 118 | REAL phis(iip1,jjp1) |
|---|
| 119 | real utime,latst,lonst,tamp,hlst |
|---|
| 120 | real localtime(iip1),heure(iip1) |
|---|
| 121 | common/temporaire/localtime |
|---|
| 122 | |
|---|
| 123 | INTEGER*4 day0 |
|---|
| 124 | c integer nmoy(jjp1,llm),tp1,tp2,pass |
|---|
| 125 | |
|---|
| 126 | |
|---|
| 127 | CHARACTER*1 yes |
|---|
| 128 | |
|---|
| 129 | LOGICAL ldrs |
|---|
| 130 | integer iformat,nid |
|---|
| 131 | |
|---|
| 132 | c declarations de l'interface avec mywrite: |
|---|
| 133 | c ----------------------------------------- |
|---|
| 134 | |
|---|
| 135 | CHARACTER file*80 |
|---|
| 136 | CHARACTER nomfich*60 |
|---|
| 137 | |
|---|
| 138 | c externe: |
|---|
| 139 | c -------- |
|---|
| 140 | |
|---|
| 141 | EXTERNAL iniconst,inigeom,covcont,mywrite |
|---|
| 142 | EXTERNAL inifilr,exner |
|---|
| 143 | EXTERNAL solarlong,coordij,moy2 |
|---|
| 144 | EXTERNAL SSUM |
|---|
| 145 | REAL SSUM |
|---|
| 146 | EXTERNAL lnblnk |
|---|
| 147 | INTEGER lnblnk |
|---|
| 148 | logical icdf,idrs |
|---|
| 149 | |
|---|
| 150 | c Dust |
|---|
| 151 | c ---- |
|---|
| 152 | |
|---|
| 153 | #include "fxyprim.h" |
|---|
| 154 | |
|---|
| 155 | c----------------------------------------------------------------------- |
|---|
| 156 | c initialisations: |
|---|
| 157 | c ---------------- |
|---|
| 158 | |
|---|
| 159 | ldrs=.true. |
|---|
| 160 | unitlec=11 |
|---|
| 161 | itautot=0 |
|---|
| 162 | idrs=.false. |
|---|
| 163 | icdf=.false. |
|---|
| 164 | |
|---|
| 165 | c Lecture du fichier a lire |
|---|
| 166 | c ------------------------- |
|---|
| 167 | |
|---|
| 168 | PRINT*,'entrer le nom du fichier NetCDF de la MCD' |
|---|
| 169 | READ(5,'(a)') nomfich |
|---|
| 170 | file=nomfich(1:lnblnk(nomfich)) |
|---|
| 171 | PRINT*,'file',file |
|---|
| 172 | |
|---|
| 173 | |
|---|
| 174 | c Lecture des autres parametres du programme |
|---|
| 175 | c ------------------------------------------ |
|---|
| 176 | |
|---|
| 177 | write(*,*) 'Raie a simuler ? ' |
|---|
| 178 | write(*,*)'[1]12CO(2-1) [2]:13CO(2-1) ', |
|---|
| 179 | & ' [3]:12CO(1-0) [4]13CO(1-0) [5]H2O(183)' |
|---|
| 180 | read(*,*) iraie |
|---|
| 181 | |
|---|
| 182 | write(*,*) 'heure local du point subterrestre ? ' |
|---|
| 183 | read(*,*) hlst |
|---|
| 184 | |
|---|
| 185 | write(*,*) 'latitude du point subterrestre ? ' |
|---|
| 186 | read(*,*) latst |
|---|
| 187 | |
|---|
| 188 | c--------------------------------------------------------- |
|---|
| 189 | c Lecture des fichiers stats , diagfi ou histmoy |
|---|
| 190 | c--------------------------------------------------------- |
|---|
| 191 | |
|---|
| 192 | 800 continue ! LOOP SUR LES FICHIERS |
|---|
| 193 | |
|---|
| 194 | c Ouverture fichiers : |
|---|
| 195 | c ------------------ |
|---|
| 196 | |
|---|
| 197 | c Ouverture fichier NetCDF |
|---|
| 198 | ierr = NF_OPEN (file(1:lnblnk(file))//'.nc', NF_NOWRITE,nid) |
|---|
| 199 | IF (ierr.NE.NF_NOERR) THEN |
|---|
| 200 | write(6,*)' Pb d''ouverture du fichier ',file |
|---|
| 201 | write(6,*)' ierr = ', ierr |
|---|
| 202 | CALL abort |
|---|
| 203 | ENDIF |
|---|
| 204 | |
|---|
| 205 | c---------------------------------------------------------------------- |
|---|
| 206 | c initialisation de la physique: |
|---|
| 207 | c ------------------------------ |
|---|
| 208 | |
|---|
| 209 | c initialisations aux valeurs martiennes |
|---|
| 210 | |
|---|
| 211 | rad=3397200. ! rayon de mars (m) ~3397200 m |
|---|
| 212 | daysec=88775. ! duree du sol (s) ~88775 s |
|---|
| 213 | omeg=4.*asin(1.)/(daysec) ! vitesse de rotation (rad.s-1) |
|---|
| 214 | g=3.72 ! gravite (m.s-2) ~3.72 |
|---|
| 215 | mugaz=43.49 ! Masse molaire de l''atm (g.mol-1) ~43.49 |
|---|
| 216 | kappa=.256793 ! = r/cp ~0.256793 |
|---|
| 217 | r = 191.18213 |
|---|
| 218 | pi=2.*asin(1.) |
|---|
| 219 | ecritphy =1 |
|---|
| 220 | iphysiq=1 |
|---|
| 221 | day_ini=0. |
|---|
| 222 | day_step=1 |
|---|
| 223 | |
|---|
| 224 | WRITE (*,*) 'day0 = ' , day0 |
|---|
| 225 | |
|---|
| 226 | CALL iniconst |
|---|
| 227 | CALL inigeom |
|---|
| 228 | CALL inifilr |
|---|
| 229 | |
|---|
| 230 | c sigma aux niveaux s: |
|---|
| 231 | |
|---|
| 232 | DO l=1,llm |
|---|
| 233 | print*,'sig_s(',l,') = ',sig_s(l) |
|---|
| 234 | ENDDO |
|---|
| 235 | |
|---|
| 236 | c Format Grads (Stats) : On suppose que l'on a 12 pas de temps : |
|---|
| 237 | c -------------------- |
|---|
| 238 | |
|---|
| 239 | nbpas=12 |
|---|
| 240 | print*,'nbpas=',nbpas |
|---|
| 241 | |
|---|
| 242 | c====================================================================== |
|---|
| 243 | c debut de la boucle sur les etats dans un fichier histoire: |
|---|
| 244 | c====================================================================== |
|---|
| 245 | |
|---|
| 246 | itau= 1 |
|---|
| 247 | |
|---|
| 248 | 801 continue ! LOOP SUR LES PAS DE TEMPS START HERE |
|---|
| 249 | |
|---|
| 250 | varname='u' |
|---|
| 251 | ierr =mcdreadnc(nid,varname,llm,u,itau) |
|---|
| 252 | varname='v' |
|---|
| 253 | ierr =mcdreadnc(nid,varname,llm,v,itau) |
|---|
| 254 | varname='temp' |
|---|
| 255 | ierr =mcdreadnc(nid,varname,llm,t,itau) |
|---|
| 256 | varname='ps' |
|---|
| 257 | ierr =mcdreadnc(nid,varname,1,ps,itau) |
|---|
| 258 | c PRINT*,'ps',ps(iip1/2,jjp1/2) |
|---|
| 259 | varname='tsurf' |
|---|
| 260 | ierr =mcdreadnc(nid,varname,1,tsurf,itau) |
|---|
| 261 | c PRINT*,'tsurf',tsurf(iip1/2,jjp1/2) |
|---|
| 262 | |
|---|
| 263 | c********************************************************** |
|---|
| 264 | c Traitement des donnees : (we are in LOOP on timestep) |
|---|
| 265 | c********************************************************** |
|---|
| 266 | |
|---|
| 267 | |
|---|
| 268 | c universal time (in stats file) |
|---|
| 269 | c------------------------------- |
|---|
| 270 | utime = (itau-1)*2. |
|---|
| 271 | |
|---|
| 272 | c Longitude du point subterrestre |
|---|
| 273 | c ------------------------------- |
|---|
| 274 | |
|---|
| 275 | tamp=(hlst-utime)*(iip1/24.) |
|---|
| 276 | if(tamp.GT.iip1) then |
|---|
| 277 | tamp=tamp-iip1 |
|---|
| 278 | endif |
|---|
| 279 | if(tamp.LT.0) then |
|---|
| 280 | tamp=iip1+tamp |
|---|
| 281 | endif |
|---|
| 282 | lonst=rlonv(nint(tamp))*180./pi |
|---|
| 283 | write(*,*) hlst,latst,lonst |
|---|
| 284 | |
|---|
| 285 | lonst=(hlst-utime)*15 |
|---|
| 286 | if(lonst.gt.180) lonst=lonst-180. |
|---|
| 287 | if(lonst.lt.-180) lonst=lonst+180. |
|---|
| 288 | write(*,*) hlst,latst,lonst |
|---|
| 289 | |
|---|
| 290 | c Local time |
|---|
| 291 | c ---------- |
|---|
| 292 | |
|---|
| 293 | c local time (in stats file) |
|---|
| 294 | DO i = 1,iim |
|---|
| 295 | localtime(i)=utime + 12.*rlonv(i)/pi |
|---|
| 296 | if(localtime(i).gt.24) localtime(i)=localtime(i)-24. |
|---|
| 297 | if(localtime(i).lt.0) localtime(i)=localtime(i)+24. |
|---|
| 298 | heure(i)=localtime(i) |
|---|
| 299 | c if(itau.EQ.1) heure(i)=0.375*(i-1) |
|---|
| 300 | c if(itau.EQ.1) print*,heure(i) |
|---|
| 301 | ENDDO |
|---|
| 302 | |
|---|
| 303 | c if(itau.EQ.1) heure(iip1)=24. |
|---|
| 304 | heure(iip1)=localtime(1) |
|---|
| 305 | |
|---|
| 306 | c ****************************************** |
|---|
| 307 | |
|---|
| 308 | |
|---|
| 309 | c c Passage en niveaux de pression |
|---|
| 310 | c -------------------------------- |
|---|
| 311 | |
|---|
| 312 | call xsig2p (ps,T,pref,llm,Tp) |
|---|
| 313 | call xsig2p (ps,u,pref,llm,up) |
|---|
| 314 | call xsig2p (ps,v,pref,llm,vp) |
|---|
| 315 | |
|---|
| 316 | c Test temporaire |
|---|
| 317 | c do l=1, llm |
|---|
| 318 | c write(77,*) pref(l),up(1,1,l) |
|---|
| 319 | c end do |
|---|
| 320 | |
|---|
| 321 | c Traitement des donnees -> creation de spectres synthetiques |
|---|
| 322 | c ----------------------------------------------------------- |
|---|
| 323 | |
|---|
| 324 | c calcul des vitesses projetees |
|---|
| 325 | |
|---|
| 326 | do l=1,llm |
|---|
| 327 | do j=1,jjp1 |
|---|
| 328 | do i=1,iip1 |
|---|
| 329 | |
|---|
| 330 | |
|---|
| 331 | up(i,j,l)=(242.*cos(rlatu(j))+up(i,j,l))* |
|---|
| 332 | . sin(rlonv(i)-lonst*pi/180)*cos(latst*pi/180) |
|---|
| 333 | |
|---|
| 334 | |
|---|
| 335 | vp(i,j,l)=vp(i,j,l)*(cos(rlonv(i)-lonst*pi/180)*sin(rlatu(j)) |
|---|
| 336 | . *cos(latst*pi/180)-sin(latst*pi/180)*cos(rlatu(j))) |
|---|
| 337 | |
|---|
| 338 | |
|---|
| 339 | enddo |
|---|
| 340 | enddo |
|---|
| 341 | enddo |
|---|
| 342 | |
|---|
| 343 | c synthese des spectres |
|---|
| 344 | |
|---|
| 345 | c print*,'itau=',itau,nbpas |
|---|
| 346 | |
|---|
| 347 | call spectre2(itau,iraie,latst,lonst,olon,olat, |
|---|
| 348 | .vp,up,tp,ps,tsurf) |
|---|
| 349 | |
|---|
| 350 | |
|---|
| 351 | |
|---|
| 352 | c Fin de la boucle sur les pas de temps ? : (we are in LOOP on timestep) |
|---|
| 353 | c **************************************** |
|---|
| 354 | |
|---|
| 355 | itau=itau+1 |
|---|
| 356 | |
|---|
| 357 | if(itau.gt.nbpas) then |
|---|
| 358 | ierr = NF_CLOSE(nid) |
|---|
| 359 | write(*,*) 'Si vous voulez entrer un nouveau fichier' |
|---|
| 360 | print*,'entrez son nom. Sinon, RETURN' |
|---|
| 361 | READ(5,'(a)') nomfich |
|---|
| 362 | if (nomfich(1:lnblnk(nomfich)).eq."") then |
|---|
| 363 | print*,'OK, on termine la' |
|---|
| 364 | else |
|---|
| 365 | file=nomfich(1:lnblnk(nomfich)) |
|---|
| 366 | PRINT*,'file',file |
|---|
| 367 | goto 800 |
|---|
| 368 | endif |
|---|
| 369 | else |
|---|
| 370 | goto 801 |
|---|
| 371 | endif |
|---|
| 372 | |
|---|
| 373 | 900 continue |
|---|
| 374 | |
|---|
| 375 | |
|---|
| 376 | 9999 PRINT*,'Fin ' |
|---|
| 377 | 1000 format(a5,3x,i4,i3,x,a39) |
|---|
| 378 | 7777 FORMAT ('latitude/longitude',4f7.1) |
|---|
| 379 | |
|---|
| 380 | END |
|---|
