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