| 1 | PROGRAM xvik |
|---|
| 2 | |
|---|
| 3 | USE filtreg_mod, ONLY: inifilr |
|---|
| 4 | USE comconst_mod, ONLY: dtvr,g,r,pi |
|---|
| 5 | |
|---|
| 6 | IMPLICIT NONE |
|---|
| 7 | c======================================================================= |
|---|
| 8 | c |
|---|
| 9 | c Pression au site Viking |
|---|
| 10 | c |
|---|
| 11 | c======================================================================= |
|---|
| 12 | c----------------------------------------------------------------------- |
|---|
| 13 | c declarations: |
|---|
| 14 | c ------------- |
|---|
| 15 | |
|---|
| 16 | |
|---|
| 17 | include "dimensions.h" |
|---|
| 18 | include "paramet.h" |
|---|
| 19 | include "comdissip.h" |
|---|
| 20 | include "comgeom2.h" |
|---|
| 21 | !#include "control.h" |
|---|
| 22 | include "netcdf.inc" |
|---|
| 23 | |
|---|
| 24 | |
|---|
| 25 | INTEGER itau,nbpas,nbpasmx |
|---|
| 26 | PARAMETER(nbpasmx=1000000) |
|---|
| 27 | REAL temps(nbpasmx) |
|---|
| 28 | INTEGER unitlec |
|---|
| 29 | INTEGER i,j,l,jj |
|---|
| 30 | REAL constR |
|---|
| 31 | |
|---|
| 32 | c Declarations NCDF: |
|---|
| 33 | c ----------------- |
|---|
| 34 | CHARACTER*100 varname |
|---|
| 35 | INTEGER ierr,nid,nvarid,dimid |
|---|
| 36 | LOGICAL nc |
|---|
| 37 | INTEGER start_ps(3),start_temp(4),start_co2ice(3) |
|---|
| 38 | INTEGER count_ps(3),count_temp(4),count_co2ice(3) |
|---|
| 39 | |
|---|
| 40 | c declarations pour les points viking: |
|---|
| 41 | c ------------------------------------ |
|---|
| 42 | INTEGER ivik(2),jvik(2),ifile(2),iv |
|---|
| 43 | REAL lonvik(2),latvik(2),phivik(2),phisim(2) |
|---|
| 44 | REAL unanj |
|---|
| 45 | |
|---|
| 46 | c variables meteo: |
|---|
| 47 | c ---------------- |
|---|
| 48 | REAL vnat(iip1,jjm,llm),unat(iip1,jjp1,llm) |
|---|
| 49 | REAL t(iip1,jjp1,llm),ps(iip1,jjp1),pstot, phis(iip1,jjp1) |
|---|
| 50 | REAL co2ice(iip1,jjp1), captotN,captotS |
|---|
| 51 | real t7(iip1,jjp1) ! temperature in 7th atmospheric layer |
|---|
| 52 | |
|---|
| 53 | REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,2),zalpha,zbeta |
|---|
| 54 | |
|---|
| 55 | LOGICAL firstcal,lcal,latcal,lvent,day_ls |
|---|
| 56 | INTEGER*4 day0 |
|---|
| 57 | |
|---|
| 58 | REAL ziceco2(iip1,jjp1) |
|---|
| 59 | REAL day,zt,sollong,sol,dayw |
|---|
| 60 | REAL airtot1,gh |
|---|
| 61 | |
|---|
| 62 | INTEGER ii,iyear,kyear |
|---|
| 63 | |
|---|
| 64 | CHARACTER*2 chr2 |
|---|
| 65 | |
|---|
| 66 | |
|---|
| 67 | c declarations de l'interface avec mywrite: |
|---|
| 68 | c ----------------------------------------- |
|---|
| 69 | |
|---|
| 70 | CHARACTER file*80 |
|---|
| 71 | CHARACTER pathchmp*80,pathsor*80,nomfich*80 |
|---|
| 72 | |
|---|
| 73 | c externe: |
|---|
| 74 | c -------- |
|---|
| 75 | |
|---|
| 76 | EXTERNAL iniconst,inigeom,covcont,mywrite |
|---|
| 77 | EXTERNAL exner,pbar |
|---|
| 78 | EXTERNAL solarlong,coordij,moy2 |
|---|
| 79 | EXTERNAL SSUM |
|---|
| 80 | REAL SSUM |
|---|
| 81 | |
|---|
| 82 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 83 | |
|---|
| 84 | c----------------------------------------------------------------------- |
|---|
| 85 | c initialisations: |
|---|
| 86 | c ---------------- |
|---|
| 87 | |
|---|
| 88 | chr2="0" |
|---|
| 89 | unanj=669. |
|---|
| 90 | print*,'WARNING!!!',unanj,'Jours/an' |
|---|
| 91 | nc=.true. |
|---|
| 92 | lcal=.true. |
|---|
| 93 | latcal=.true. |
|---|
| 94 | lvent=.false. |
|---|
| 95 | day_ls=.true. |
|---|
| 96 | |
|---|
| 97 | c lecture du fichier xvik.def |
|---|
| 98 | |
|---|
| 99 | phivik(1)=-3637 |
|---|
| 100 | phivik(2)=-4505 |
|---|
| 101 | |
|---|
| 102 | |
|---|
| 103 | |
|---|
| 104 | OPEN(99,file='xvik.def',form='formatted') |
|---|
| 105 | |
|---|
| 106 | READ(99,*) |
|---|
| 107 | READ(99,*,iostat=ierr) phivik |
|---|
| 108 | IF(ierr.NE.0) GOTO 105 |
|---|
| 109 | |
|---|
| 110 | READ(99,*,END=105) |
|---|
| 111 | READ(99,'(a)',END=105) pathchmp |
|---|
| 112 | READ(99,*,END=105) |
|---|
| 113 | READ(99,'(a)',END=105) pathsor |
|---|
| 114 | READ(99,*,END=105) |
|---|
| 115 | c READ(99,'(l1)',END=105) day_ls |
|---|
| 116 | READ(99,'(l1)',END=105) |
|---|
| 117 | READ(99,'(l1)',END=105) lcal |
|---|
| 118 | READ(99,'(l1)',END=105) |
|---|
| 119 | READ(99,'(l1)',END=105) lvent |
|---|
| 120 | READ(99,'(l1)',END=105) |
|---|
| 121 | READ(99,'(l1)',END=105) latcal |
|---|
| 122 | |
|---|
| 123 | 105 CONTINUE |
|---|
| 124 | CLOSE(99) |
|---|
| 125 | write (*,*)'>>>>>>>>>>>>>>>>', phivik,g |
|---|
| 126 | DO iv=1,2 |
|---|
| 127 | phivik(iv)=phivik(iv)*3.73 |
|---|
| 128 | END DO |
|---|
| 129 | |
|---|
| 130 | write(*,*) ' pathchmp:',trim(pathchmp) |
|---|
| 131 | write(*,*) ' pathsor:',trim(pathsor) |
|---|
| 132 | |
|---|
| 133 | c----------------------------------------------------------------------- |
|---|
| 134 | c----------------------------------------------------------------------- |
|---|
| 135 | c ouverture des fichiers xgraph: |
|---|
| 136 | c ------------------------------ |
|---|
| 137 | |
|---|
| 138 | ifile(1)=12 |
|---|
| 139 | ifile(2)=13 |
|---|
| 140 | kyear=-1 |
|---|
| 141 | c OPEN(77,file='xlongday',form='formatted') |
|---|
| 142 | |
|---|
| 143 | unitlec=11 |
|---|
| 144 | c |
|---|
| 145 | PRINT*,'entrer le nom du fichier NC' |
|---|
| 146 | READ(5,'(a)') nomfich |
|---|
| 147 | |
|---|
| 148 | PRINT *,'nomfich : ',nomfich |
|---|
| 149 | |
|---|
| 150 | |
|---|
| 151 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 152 | c grande boucle sur les fichiers histoire: |
|---|
| 153 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 154 | |
|---|
| 155 | firstcal=.true. |
|---|
| 156 | DO WHILE(len_trim(nomfich).GT.0.AND.len_trim(nomfich).LT.50) |
|---|
| 157 | PRINT *,'>>> nomfich : ',trim(nomfich) |
|---|
| 158 | |
|---|
| 159 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 160 | |
|---|
| 161 | file=pathchmp(1:len_trim(pathchmp))//'/'// |
|---|
| 162 | s nomfich(1:len_trim(nomfich)) |
|---|
| 163 | PRINT*,'file.nc: ', file(1:len_trim(file))//'.nc' |
|---|
| 164 | PRINT*,'timestep ',dtvr |
|---|
| 165 | |
|---|
| 166 | IF(nc) THEN |
|---|
| 167 | ierr= NF_OPEN(file(1:len_trim(file))//'.nc',NF_NOWRITE,nid) |
|---|
| 168 | ELSE |
|---|
| 169 | PRINT*,'Ouverture binaire ',file |
|---|
| 170 | OPEN(unitlec,file=file,status='old',form='unformatted', |
|---|
| 171 | . iostat=ierr) |
|---|
| 172 | ENDIF |
|---|
| 173 | |
|---|
| 174 | c---------------------------------------------------------------------- |
|---|
| 175 | c initialisation de la physique: |
|---|
| 176 | c ------------------------------ |
|---|
| 177 | |
|---|
| 178 | CALL readhead_NC(file(1:len_trim(file))//'.nc',day0,phis,constR) |
|---|
| 179 | |
|---|
| 180 | WRITE (*,*) 'day0 = ' , day0 |
|---|
| 181 | |
|---|
| 182 | CALL conf_gcm( 99, .TRUE. ) |
|---|
| 183 | CALL iniconst |
|---|
| 184 | CALL inigeom |
|---|
| 185 | ! CALL inifilr |
|---|
| 186 | |
|---|
| 187 | |
|---|
| 188 | c Lecture temps : |
|---|
| 189 | |
|---|
| 190 | ierr= NF_INQ_DIMID (nid,"Time",dimid) |
|---|
| 191 | IF (ierr.NE.NF_NOERR) THEN |
|---|
| 192 | PRINT*, 'xvik: Le champ <Time> est absent' |
|---|
| 193 | CALL abort |
|---|
| 194 | ENDIF |
|---|
| 195 | |
|---|
| 196 | ierr= NF_INQ_DIMLEN (nid,dimid,nbpas) |
|---|
| 197 | |
|---|
| 198 | ierr = NF_INQ_VARID (nid, "Time", nvarid) |
|---|
| 199 | #ifdef NC_DOUBLE |
|---|
| 200 | ierr = NF_GET_VAR_DOUBLE(nid, nvarid, temps) |
|---|
| 201 | #else |
|---|
| 202 | ierr = NF_GET_VAR_REAL(nid, nvarid, temps) |
|---|
| 203 | #endif |
|---|
| 204 | IF (ierr.NE.NF_NOERR) THEN |
|---|
| 205 | PRINT*, 'xvik: Lecture echouee pour <Time>' |
|---|
| 206 | CALL abort |
|---|
| 207 | ENDIF |
|---|
| 208 | |
|---|
| 209 | PRINT*,'temps(1:10)',(temps(itau),itau=1,10) |
|---|
| 210 | |
|---|
| 211 | c----------------------------------------------------------------------- |
|---|
| 212 | c coordonnees des point Viking: |
|---|
| 213 | c ----------------------------- |
|---|
| 214 | |
|---|
| 215 | latvik(1)=22.27*pi/180. |
|---|
| 216 | lonvik(1)=-47.95*pi/180. |
|---|
| 217 | latvik(2)=47.67*pi/180. |
|---|
| 218 | lonvik(2)=(360.-225.71)*pi/180. |
|---|
| 219 | |
|---|
| 220 | c ponderations pour les 4 points autour de Viking |
|---|
| 221 | DO iv=1,2 |
|---|
| 222 | ! locate index of GCM grid points near VL |
|---|
| 223 | do i=1,iim |
|---|
| 224 | ! we know longitudes are ordered -180...180 |
|---|
| 225 | if ((lonvik(iv).ge.rlonu(i)).and. |
|---|
| 226 | & (lonvik(iv).le.rlonu(i+1))) then |
|---|
| 227 | ivik(iv)=i |
|---|
| 228 | exit |
|---|
| 229 | endif |
|---|
| 230 | enddo |
|---|
| 231 | do j=1,jjm-1 |
|---|
| 232 | !we know tha latitudes are ordered 90...-90 |
|---|
| 233 | if ((latvik(iv).le.rlatv(j)).and. |
|---|
| 234 | & (latvik(iv).ge.rlatv(j+1))) then |
|---|
| 235 | jvik(iv)=j |
|---|
| 236 | exit |
|---|
| 237 | endif |
|---|
| 238 | enddo |
|---|
| 239 | zalpha=(lonvik(iv)-rlonu(ivik(iv)))/ |
|---|
| 240 | s (rlonu(ivik(iv)+1)-rlonu(ivik(iv))) |
|---|
| 241 | zbeta=(latvik(iv)-rlatv(jvik(iv)))/ |
|---|
| 242 | s (rlatv(jvik(iv)+1)-rlatv(jvik(iv))) |
|---|
| 243 | zw(0,0,iv)=(1.-zalpha)*(1.-zbeta) |
|---|
| 244 | zw(1,0,iv)=zalpha*(1.-zbeta) |
|---|
| 245 | zw(0,1,iv)=(1.-zalpha)*zbeta |
|---|
| 246 | zw(1,1,iv)=zalpha*zbeta |
|---|
| 247 | ENDDO |
|---|
| 248 | |
|---|
| 249 | c altitude reelle et modele aux points Viking |
|---|
| 250 | DO iv=1,2 |
|---|
| 251 | phisim(iv)=0. |
|---|
| 252 | DO jj=0,1 |
|---|
| 253 | j=jvik(iv)+jj |
|---|
| 254 | DO ii=0,1 |
|---|
| 255 | i=ivik(iv)+ii |
|---|
| 256 | phisim(iv)=phisim(iv)+zw(ii,jj,iv)*phis(i,j) |
|---|
| 257 | ENDDO |
|---|
| 258 | ENDDO |
|---|
| 259 | ENDDO |
|---|
| 260 | PRINT*,'relief aux points Viking pour les sorties:',phivik |
|---|
| 261 | |
|---|
| 262 | c---------------------------------------------------------------------- |
|---|
| 263 | c lectures des etats: |
|---|
| 264 | c ------------------- |
|---|
| 265 | |
|---|
| 266 | airtot1=1./(SSUM(ip1jmp1,aire,1)-SSUM(jjp1,aire,iip1)) |
|---|
| 267 | |
|---|
| 268 | c====================================================================== |
|---|
| 269 | c debut de la boucle sur les etats dans un fichier histoire: |
|---|
| 270 | c====================================================================== |
|---|
| 271 | count_ps=(/iip1,jjp1,1/) |
|---|
| 272 | count_co2ice=(/iip1,jjp1,1/) |
|---|
| 273 | count_temp=(/iip1,jjp1,llm,1/) |
|---|
| 274 | |
|---|
| 275 | DO itau=1,nbpas |
|---|
| 276 | |
|---|
| 277 | start_ps=(/1,1,itau/) |
|---|
| 278 | start_co2ice=(/1,1,itau/) |
|---|
| 279 | start_temp=(/1,1,1,itau/) |
|---|
| 280 | c lecture drs des champs: |
|---|
| 281 | c ----------------------- |
|---|
| 282 | c varname='u' |
|---|
| 283 | c ierr=drsread (unitlec,varname,unat,itau) |
|---|
| 284 | c PRINT*,'unat',unat(iip1/2,jjp1/2,llm/2) |
|---|
| 285 | c varname='v' |
|---|
| 286 | c ierr=drsread (unitlec,varname,vnat,itau) |
|---|
| 287 | c PRINT*,'vnat',vnat(iip1/2,jjp1/2,llm/2) |
|---|
| 288 | |
|---|
| 289 | ccccccccc LECTURE Ps ccccccccccccccccccccccccccc |
|---|
| 290 | ierr = NF_INQ_VARID (nid, "ps", nvarid) |
|---|
| 291 | #ifdef NC_DOUBLE |
|---|
| 292 | ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start_ps,count_ps, ps) |
|---|
| 293 | #else |
|---|
| 294 | ierr = NF_GET_VARA_REAL(nid, nvarid,start_ps,count_ps, ps) |
|---|
| 295 | #endif |
|---|
| 296 | IF (ierr.NE.NF_NOERR) THEN |
|---|
| 297 | PRINT*, 'xvik: Lecture echouee pour <ps>' |
|---|
| 298 | CALL abort |
|---|
| 299 | ENDIF |
|---|
| 300 | |
|---|
| 301 | PRINT*,'ps',ps(iip1/2,jjp1/2) |
|---|
| 302 | |
|---|
| 303 | ccccccccc LECTURE Temperature ccccccccccccccccccccccccccc |
|---|
| 304 | ierr = NF_INQ_VARID (nid, "temp", nvarid) |
|---|
| 305 | #ifdef NC_DOUBLE |
|---|
| 306 | ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_temp,count_temp, t) |
|---|
| 307 | #else |
|---|
| 308 | ierr = NF_GET_VARA_REAL(nid,nvarid,start_temp,count_temp, t) |
|---|
| 309 | #endif |
|---|
| 310 | IF (ierr.NE.NF_NOERR) THEN |
|---|
| 311 | PRINT*, 'xvik: Lecture echouee pour <temp>' |
|---|
| 312 | ! Ehouarn: proceed anyways |
|---|
| 313 | ! CALL abort |
|---|
| 314 | write(*,*)'--> Setting temperature to zero !!!' |
|---|
| 315 | t(1:iip1,1:jjp1,1:llm)=0.0 |
|---|
| 316 | write(*,*)'--> looking for temp7 (temp in 7th layer)' |
|---|
| 317 | ierr=NF_INQ_VARID(nid,"temp7", nvarid) |
|---|
| 318 | if (ierr.eq.NF_NOERR) then |
|---|
| 319 | write(*,*) " OK, found temp7 variable" |
|---|
| 320 | #ifdef NC_DOUBLE |
|---|
| 321 | ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start_ps,count_ps,t7) |
|---|
| 322 | #else |
|---|
| 323 | ierr=NF_GET_VARA_REAL(nid,nvarid,start_ps,count_ps,t7) |
|---|
| 324 | #endif |
|---|
| 325 | if (ierr.ne.NF_NOERR) then |
|---|
| 326 | write(*,*)'xvik: failed loading temp7 !' |
|---|
| 327 | stop |
|---|
| 328 | endif |
|---|
| 329 | else ! no 'temp7' variable |
|---|
| 330 | write(*,*)' No temp7 variable either !' |
|---|
| 331 | write(*,*)' Will have to to without ...' |
|---|
| 332 | t7(1:iip1,1:jjp1)=0.0 |
|---|
| 333 | endif |
|---|
| 334 | ELSE ! t() was successfully loaded, copy 7th layer to t7() |
|---|
| 335 | t7(1:iip1,1:jjp1)=t(1:iip1,1:jjp1,7) |
|---|
| 336 | ENDIF |
|---|
| 337 | |
|---|
| 338 | c PRINT*,'t',t(iip1/2,jjp1/2,llm/2) |
|---|
| 339 | |
|---|
| 340 | ccccccccc LECTURE co2ice ccccccccccccccccccccccccccc |
|---|
| 341 | ierr = NF_INQ_VARID (nid, "co2ice", nvarid) |
|---|
| 342 | #ifdef NC_DOUBLE |
|---|
| 343 | ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_co2ice, |
|---|
| 344 | & count_co2ice, co2ice) |
|---|
| 345 | #else |
|---|
| 346 | ierr = NF_GET_VARA_REAL(nid, nvarid,start_co2ice, |
|---|
| 347 | & count_co2ice, co2ice) |
|---|
| 348 | #endif |
|---|
| 349 | IF (ierr.NE.NF_NOERR) THEN |
|---|
| 350 | PRINT*, 'xvik: Lecture echouee pour <co2ice>' |
|---|
| 351 | CALL abort |
|---|
| 352 | ENDIF |
|---|
| 353 | |
|---|
| 354 | |
|---|
| 355 | c Gestion du temps |
|---|
| 356 | c ---------------- |
|---|
| 357 | day=temps(itau) |
|---|
| 358 | PRINT*,'day ',day |
|---|
| 359 | CALL solarlong(day+day0,sollong) |
|---|
| 360 | sol=day+day0+461. ! aded offset to sync with VL mission "sol 1" |
|---|
| 361 | iyear=sol/unanj |
|---|
| 362 | WRITE (*,*) 'iyear',iyear |
|---|
| 363 | sol=sol-iyear*unanj |
|---|
| 364 | c |
|---|
| 365 | c Ouverture / fermeture des fichiers |
|---|
| 366 | c ---------------------------------- |
|---|
| 367 | IF (iyear.NE.kyear) THEN |
|---|
| 368 | WRITE(chr2(1:1),'(i1)') iyear+1 |
|---|
| 369 | WRITE (*,*) 'iyear bis',iyear |
|---|
| 370 | WRITE (*,*) 'chr2' |
|---|
| 371 | WRITE (*,*) chr2 |
|---|
| 372 | IF(iyear.GE.9) WRITE(chr2,'(i2)') iyear+1 |
|---|
| 373 | kyear=iyear |
|---|
| 374 | DO ii=1,2 |
|---|
| 375 | CLOSE(10+ifile(ii)) |
|---|
| 376 | CLOSE(2+ifile(ii)) |
|---|
| 377 | CLOSE(4+ifile(ii)) |
|---|
| 378 | CLOSE(6+ifile(ii)) |
|---|
| 379 | CLOSE(8+ifile(ii)) |
|---|
| 380 | CLOSE(16+ifile(ii)) |
|---|
| 381 | CLOSE(12+ifile(ii)) |
|---|
| 382 | CLOSE(14+ifile(ii)) |
|---|
| 383 | CLOSE(97) |
|---|
| 384 | CLOSE(98) |
|---|
| 385 | ENDDO |
|---|
| 386 | CLOSE(5+ifile(1)) |
|---|
| 387 | OPEN(ifile(1)+10,file='xpsol1'//chr2,form='formatted') |
|---|
| 388 | OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted') |
|---|
| 389 | c OPEN(ifile(1)+8,file='xbpsol1'//chr2,form='formatted') |
|---|
| 390 | c OPEN(ifile(2)+8,file='xbpsol2'//chr2,form='formatted') |
|---|
| 391 | c OPEN(ifile(1)+2,file='xlps1'//chr2,form='formatted') |
|---|
| 392 | c OPEN(ifile(2)+2,file='xlps2'//chr2,form='formatted') |
|---|
| 393 | IF(lcal) THEN |
|---|
| 394 | c OPEN(ifile(2)+4,file='xpressud'//chr2,form='formatted') |
|---|
| 395 | c OPEN(ifile(1)+4,file='xpresnord'//chr2,form='formatted') |
|---|
| 396 | c OPEN(ifile(1)+6,file='xpm2'//chr2,form='formatted') |
|---|
| 397 | ENDIF |
|---|
| 398 | IF(latcal) THEN |
|---|
| 399 | c OPEN(ifile(2)+14,file='xlats'//chr2,form='formatted') |
|---|
| 400 | c OPEN(ifile(1)+14,file='xlatn'//chr2,form='formatted') |
|---|
| 401 | ENDIF |
|---|
| 402 | IF(lvent) THEN |
|---|
| 403 | c OPEN(ifile(1)+16,file='xu1'//chr2,form='formatted') |
|---|
| 404 | c OPEN(ifile(2)+16,file='xu2'//chr2,form='formatted') |
|---|
| 405 | c OPEN(ifile(1)+12,file='xv1'//chr2,form='formatted') |
|---|
| 406 | c OPEN(ifile(2)+12,file='xv2'//chr2,form='formatted') |
|---|
| 407 | ENDIF |
|---|
| 408 | OPEN(97,file='xprestot'//chr2,form='formatted') |
|---|
| 409 | c OPEN(98,file='xlat37_'//chr2,form='formatted') |
|---|
| 410 | c WRITE(98,'(f5.1,16f7.1)') 0.,(rlonv(i)*180./pi,i=1,iim,4) |
|---|
| 411 | ENDIF |
|---|
| 412 | |
|---|
| 413 | |
|---|
| 414 | sollong=sollong*180./pi |
|---|
| 415 | IF(day_ls) THEN |
|---|
| 416 | dayw=sol |
|---|
| 417 | write(*,*) 'dayw', dayw |
|---|
| 418 | ELSE |
|---|
| 419 | dayw=sollong |
|---|
| 420 | ENDIF |
|---|
| 421 | |
|---|
| 422 | c Calcul de la moyenne planetaire |
|---|
| 423 | c ------------------------------- |
|---|
| 424 | pstot=0. |
|---|
| 425 | captotS=0. |
|---|
| 426 | captotN=0. |
|---|
| 427 | DO j=1,jjp1 |
|---|
| 428 | DO i=1,iim |
|---|
| 429 | pstot=pstot+aire(i,j)*ps(i,j) |
|---|
| 430 | ENDDO |
|---|
| 431 | ENDDO |
|---|
| 432 | |
|---|
| 433 | DO j=1,jjp1/2 |
|---|
| 434 | DO i=1,iim |
|---|
| 435 | captotN = captotN +aire(i,j)*co2ice(i,j) |
|---|
| 436 | ENDDO |
|---|
| 437 | ENDDO |
|---|
| 438 | DO j=jjp1/2+1, jjp1 |
|---|
| 439 | DO i=1,iim |
|---|
| 440 | captotS = captotS +aire(i,j)*co2ice(i,j) |
|---|
| 441 | ENDDO |
|---|
| 442 | ENDDO |
|---|
| 443 | WRITE(97,'(4e16.6)') dayw,pstot*airtot1 |
|---|
| 444 | & , captotN*g*airtot1, captotS*g*airtot1 |
|---|
| 445 | |
|---|
| 446 | IF(.NOT.firstcal) THEN |
|---|
| 447 | c WRITE(98,'(f5.1,16f7.3)') |
|---|
| 448 | c s dayw,(ps(i,37),i=1,iim,4) |
|---|
| 449 | |
|---|
| 450 | c boucle sur les sites vikings: |
|---|
| 451 | c ---------------------------- |
|---|
| 452 | |
|---|
| 453 | DO iv=1,2 |
|---|
| 454 | |
|---|
| 455 | c interpolation de la temperature dans la 7eme couche, de la pression |
|---|
| 456 | c de surface et des vents aux points viking. |
|---|
| 457 | |
|---|
| 458 | zp1=0. |
|---|
| 459 | zp2=0. |
|---|
| 460 | zp2_sm=0. |
|---|
| 461 | zt=0. |
|---|
| 462 | ! zu=0. |
|---|
| 463 | ! zv=0. |
|---|
| 464 | DO jj=0,1 |
|---|
| 465 | j=jvik(iv)+jj |
|---|
| 466 | DO ii=0,1 |
|---|
| 467 | i=ivik(iv)+ii |
|---|
| 468 | ! zt=zt+zw(ii,jj,iv)*t(i,j,7) |
|---|
| 469 | zt=zt+zw(ii,jj,iv)*t7(i,j) |
|---|
| 470 | ! zp1=zp1+zw(ii,jj,iv)*ps(i,j) |
|---|
| 471 | zp1=zp1+zw(ii,jj,iv)*log(ps(i,j)) ! interpolate in log(P) |
|---|
| 472 | WRITE (*,*) 'ps autour iv',ps(i,j),iv |
|---|
| 473 | ! zu=zu+zw(ii,jj,iv)*unat(i,j,1)/cu(i,j) |
|---|
| 474 | ! zv=zv+zw(ii,jj,iv)*vnat(i,j,1)/cv(i,j) |
|---|
| 475 | ENDDO |
|---|
| 476 | ENDDO |
|---|
| 477 | zp1=exp(zp1) ! because of the bilinear interpolation in log(P) |
|---|
| 478 | |
|---|
| 479 | c pression au sol extrapolee a partir de la temp. 7eme couche |
|---|
| 480 | WRITE (*,*) 'constR ',constR |
|---|
| 481 | WRITE (*,*) 'zt ',zt |
|---|
| 482 | gh=constR*zt |
|---|
| 483 | if (gh.eq.0) then ! if we don't have temperature values |
|---|
| 484 | ! assume a scale height of 10km |
|---|
| 485 | zp2=zp1*exp(-(phivik(iv)-phisim(iv))/(3.73*1.e4)) |
|---|
| 486 | else |
|---|
| 487 | zp2=zp1*exp(-(phivik(iv)-phisim(iv))/gh) |
|---|
| 488 | endif |
|---|
| 489 | WRITE (*,*) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh' |
|---|
| 490 | WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh |
|---|
| 491 | ! WRITE(ifile(iv)+10,'(2e15.5)') dayw,zp1 |
|---|
| 492 | WRITE(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1 |
|---|
| 493 | |
|---|
| 494 | c sorties eventuelles de vent |
|---|
| 495 | ! IF(lvent) THEN |
|---|
| 496 | ! WRITE(ifile(iv)+16,'(2e15.5)') |
|---|
| 497 | ! s dayw,zu |
|---|
| 498 | ! WRITE(ifile(iv)+12,'(2e15.5)') |
|---|
| 499 | ! s dayw,zv |
|---|
| 500 | ! ENDIF |
|---|
| 501 | ENDDO |
|---|
| 502 | c IF (lcal) THEN |
|---|
| 503 | c WRITE(ifile(1)+4,'(2e15.6)') dayw,airtot1*g*.01* |
|---|
| 504 | c s (SSUM(ip1jmp1/2,ziceco2,1)-SSUM(jjp1/2,ziceco2,iip1)) |
|---|
| 505 | c WRITE(ifile(2)+4,'(2e15.6)') dayw,airtot1*g*.01* |
|---|
| 506 | c s (SSUM(iip1*jjm/2,ziceco2(1,jjm/2+2),1)- |
|---|
| 507 | c s SSUM(jjm/2,ziceco2(1,jjm/2+2),iip1)) |
|---|
| 508 | c ENDIF |
|---|
| 509 | c IF(latcal) THEN |
|---|
| 510 | c CALL icelat(iim,jjm,ziceco2,rlatv,zicelat) |
|---|
| 511 | c WRITE(ifile(1)+14,'(2e15.6)') dayw,zicelat(1)*180./pi |
|---|
| 512 | c WRITE(ifile(2)+14,'(2e15.6)') dayw,zicelat(2)*180./pi |
|---|
| 513 | c ENDIF |
|---|
| 514 | ENDIF |
|---|
| 515 | firstcal=.false. |
|---|
| 516 | |
|---|
| 517 | c====================================================================== |
|---|
| 518 | c Fin de la boucle sur les etats du fichier histoire: |
|---|
| 519 | c====================================================================== |
|---|
| 520 | ENDDO |
|---|
| 521 | |
|---|
| 522 | ierr= NF_CLOSE(nid) |
|---|
| 523 | |
|---|
| 524 | PRINT*,'Fin du fichier',nomfich |
|---|
| 525 | print*,'Entrer un nouveau fichier ou return pour finir' |
|---|
| 526 | READ(5,'(a)',err=9999) nomfich |
|---|
| 527 | |
|---|
| 528 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 529 | c Fin de la boucle sur les fichiers histoire: |
|---|
| 530 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 531 | |
|---|
| 532 | ENDDO |
|---|
| 533 | |
|---|
| 534 | PRINT*,'relief du point V1',.001*phis(ivik(1),jvik(1))/g |
|---|
| 535 | PRINT*,'relief du point V2',.001*phis(ivik(2),jvik(2))/g |
|---|
| 536 | DO iv=1,2 |
|---|
| 537 | PRINT*,'Viking',iv,' i=',ivik(iv),'j =',jvik(iv) |
|---|
| 538 | WRITE(6,7777) |
|---|
| 539 | s (rlonv(i)*180./pi,i=ivik(iv)-1,ivik(iv)+2) |
|---|
| 540 | print* |
|---|
| 541 | DO j=jvik(iv)-1,jvik(iv)+2 |
|---|
| 542 | WRITE(6,'(f8.1,10x,5f7.1)') |
|---|
| 543 | s rlatu(j)*180./pi,(phis(i,j)/(g*1000.),i=ivik(iv)-1,ivik(iv)+2) |
|---|
| 544 | ENDDO |
|---|
| 545 | print* |
|---|
| 546 | print*,'zw' |
|---|
| 547 | write(6,'(2(2f10.4/))') ((zw(ii,jj,iv),ii=0,1),jj=0,1) |
|---|
| 548 | print*,'altitude interpolee (km) ',phisim(iv)/1000./g |
|---|
| 549 | ENDDO |
|---|
| 550 | PRINT*,'R=',r |
|---|
| 551 | 9999 PRINT*,'Fin ' |
|---|
| 552 | |
|---|
| 553 | 7777 FORMAT ('latitude/longitude',4f7.1) |
|---|
| 554 | END |
|---|
