| 1 | c======================================================================= |
|---|
| 2 | SUBROUTINE datareadnc(relief,phisinit,alb,ith,z0, |
|---|
| 3 | & zmea,zstd,zsig,zgam,zthe, |
|---|
| 4 | & hmons,summit,base,zavg) |
|---|
| 5 | c======================================================================= |
|---|
| 6 | c |
|---|
| 7 | c |
|---|
| 8 | c Author: F. Hourdin 01/1997 |
|---|
| 9 | c ------- |
|---|
| 10 | c |
|---|
| 11 | c Object: To read data from Martian surface to use in a GCM |
|---|
| 12 | c ------ from NetCDF file "surface.nc" |
|---|
| 13 | c |
|---|
| 14 | c |
|---|
| 15 | c Arguments: |
|---|
| 16 | c ---------- |
|---|
| 17 | c |
|---|
| 18 | c Inputs: |
|---|
| 19 | c ------ |
|---|
| 20 | c |
|---|
| 21 | c Outputs: |
|---|
| 22 | c -------- |
|---|
| 23 | c |
|---|
| 24 | c======================================================================= |
|---|
| 25 | c donnees ALBEDO, INERTIE THERMIQUE, RELIEF: |
|---|
| 26 | c |
|---|
| 27 | c Ces donnees sont au format NetCDF dans le fichier "surface.nc" |
|---|
| 28 | c |
|---|
| 29 | c 360 valeurs en longitude (de -179.5 a 179.5) |
|---|
| 30 | c 180 valeurs en latitudes (de 89.5 a -89.5) |
|---|
| 31 | c |
|---|
| 32 | c Pour les passer au format de la grille, on utilise "interp_horiz.F" |
|---|
| 33 | c |
|---|
| 34 | c Il faut donc que ces donnees soient au format grille scalaire |
|---|
| 35 | c (imold+1 jmold+1) |
|---|
| 36 | c avec passage des coordonnees de la "boite" (rlonu, rlatv) |
|---|
| 37 | c |
|---|
| 38 | c On prend imd (d pour donnees!) |
|---|
| 39 | c imd = 360 avec copie de la 1ere valeur sur la imd+1 |
|---|
| 40 | c (rlonud de -179 a -181) |
|---|
| 41 | c jmd = 179 |
|---|
| 42 | c (rlatvd de 89 a -89) |
|---|
| 43 | c======================================================================= |
|---|
| 44 | |
|---|
| 45 | ! to use 'getin' |
|---|
| 46 | use ioipsl_getincom, only: getin |
|---|
| 47 | use comconst_mod, only: g,pi |
|---|
| 48 | use datafile_mod, only: datadir |
|---|
| 49 | use avg_horiz_mod, only: avg_horiz |
|---|
| 50 | use mvc_horiz_mod, only: mvc_horiz |
|---|
| 51 | |
|---|
| 52 | implicit none |
|---|
| 53 | |
|---|
| 54 | include "dimensions.h" |
|---|
| 55 | include "paramet.h" |
|---|
| 56 | include "comgeom.h" |
|---|
| 57 | include "netcdf.inc" |
|---|
| 58 | |
|---|
| 59 | c======================================================================= |
|---|
| 60 | c Declarations: |
|---|
| 61 | C======================================================================= |
|---|
| 62 | |
|---|
| 63 | INTEGER imd,jmd,imdp1,jmdp1 |
|---|
| 64 | parameter (imd=360,jmd=179,imdp1=361,jmdp1=180) |
|---|
| 65 | |
|---|
| 66 | INTEGER iimp1 |
|---|
| 67 | parameter (iimp1=iim+1-1/iim) |
|---|
| 68 | |
|---|
| 69 | ! Arguments: |
|---|
| 70 | CHARACTER(len=3),intent(inout) :: relief |
|---|
| 71 | REAL,intent(out) :: phisinit(iimp1*jjp1) |
|---|
| 72 | REAL,intent(out) :: alb(iimp1*jjp1) |
|---|
| 73 | REAL,intent(out) :: ith(iimp1*jjp1) |
|---|
| 74 | REAL,intent(out) :: z0(iimp1*jjp1) |
|---|
| 75 | REAL,intent(out) :: zmea(imdp1*jmdp1) |
|---|
| 76 | REAL,intent(out) :: zstd(imdp1*jmdp1) |
|---|
| 77 | REAL,intent(out) :: zsig(imdp1*jmdp1) |
|---|
| 78 | REAL,intent(out) :: zgam(imdp1*jmdp1) |
|---|
| 79 | REAL,intent(out) :: zthe(imdp1*jmdp1) |
|---|
| 80 | REAL,intent(out) :: hmons(imdp1*jmdp1) !CW17,361*180 hmons |
|---|
| 81 | REAL,intent(out) :: summit(imdp1*jmdp1) |
|---|
| 82 | REAL,intent(out) :: base(imdp1*jmdp1) |
|---|
| 83 | REAL,intent(out) :: zavg(imdp1*jmdp1) |
|---|
| 84 | |
|---|
| 85 | ! Local variables: |
|---|
| 86 | REAL zdata(imd*jmdp1) |
|---|
| 87 | REAL zdataS(imdp1*jmdp1) |
|---|
| 88 | REAL pfield(iimp1*jjp1) |
|---|
| 89 | |
|---|
| 90 | INTEGER ierr |
|---|
| 91 | |
|---|
| 92 | INTEGER unit,nvarid |
|---|
| 93 | |
|---|
| 94 | INTEGER i,j,k |
|---|
| 95 | |
|---|
| 96 | INTEGER klatdat,ngridmxgdat |
|---|
| 97 | PARAMETER (klatdat=180,ngridmxgdat=360) |
|---|
| 98 | |
|---|
| 99 | c on passe une grille en rlonu rlatv et im+1 jm a interp_horiz) |
|---|
| 100 | |
|---|
| 101 | REAL longitude(imd),latitude(jmdp1) ! Pour lecture des donnees |
|---|
| 102 | REAL rlonud(imdp1),rlatvd(jmd) |
|---|
| 103 | |
|---|
| 104 | CHARACTER*20 string |
|---|
| 105 | DIMENSION string(0:7) |
|---|
| 106 | |
|---|
| 107 | |
|---|
| 108 | !#include "lmdstd.h" |
|---|
| 109 | !#include "fxyprim.h" |
|---|
| 110 | |
|---|
| 111 | pi=2.*ASIN(1.) |
|---|
| 112 | |
|---|
| 113 | c======================================================================= |
|---|
| 114 | c rlonud, rlatvd |
|---|
| 115 | c======================================================================= |
|---|
| 116 | |
|---|
| 117 | c----------------------------------------------------------------------- |
|---|
| 118 | c Lecture NetCDF des donnees latitude et longitude |
|---|
| 119 | c----------------------------------------------------------------------- |
|---|
| 120 | write(*,*) 'datareadnc: opening file surface.nc' |
|---|
| 121 | |
|---|
| 122 | datadir="/u/lmdz/WWW/planets/mars/datadir" ! default path to surface.nc |
|---|
| 123 | call getin("datadir",datadir) ! but users may specify another path |
|---|
| 124 | |
|---|
| 125 | ierr = NF_OPEN (trim(datadir)//'/surface.nc', |
|---|
| 126 | & NF_NOWRITE,unit) |
|---|
| 127 | IF (ierr.NE.NF_NOERR) THEN |
|---|
| 128 | write(*,*)'Error : cannot open file surface.nc ' |
|---|
| 129 | write(*,*)'(in phymars/datareadnc.F)' |
|---|
| 130 | write(*,*)'It should be in :',trim(datadir),'/' |
|---|
| 131 | write(*,*)'1) You can set this path in the |
|---|
| 132 | & callphys.def file:' |
|---|
| 133 | write(*,*)' datadir=/path/to/the/datafiles' |
|---|
| 134 | write(*,*)'2) If necessary, surface.nc (and other datafiles)' |
|---|
| 135 | write(*,*)' can be obtained online on:' |
|---|
| 136 | write(*,*)'http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir' |
|---|
| 137 | CALL ABORT |
|---|
| 138 | ENDIF |
|---|
| 139 | |
|---|
| 140 | c |
|---|
| 141 | c Lecture des latitudes (coordonnees): |
|---|
| 142 | c |
|---|
| 143 | ierr = NF_INQ_VARID (unit, "latitude", nvarid) |
|---|
| 144 | #ifdef NC_DOUBLE |
|---|
| 145 | ierr = NF_GET_VAR_DOUBLE(unit, nvarid, latitude) |
|---|
| 146 | #else |
|---|
| 147 | ierr = NF_GET_VAR_REAL(unit, nvarid, latitude) |
|---|
| 148 | #endif |
|---|
| 149 | c |
|---|
| 150 | c Lecture des longitudes (coordonnees): |
|---|
| 151 | c |
|---|
| 152 | ierr = NF_INQ_VARID (unit, "longitude", nvarid) |
|---|
| 153 | #ifdef NC_DOUBLE |
|---|
| 154 | ierr = NF_GET_VAR_DOUBLE(unit, nvarid, longitude) |
|---|
| 155 | #else |
|---|
| 156 | ierr = NF_GET_VAR_REAL(unit, nvarid, longitude) |
|---|
| 157 | #endif |
|---|
| 158 | |
|---|
| 159 | c----------------------------------------------------------------------- |
|---|
| 160 | c Passage au format boites scalaires |
|---|
| 161 | c----------------------------------------------------------------------- |
|---|
| 162 | |
|---|
| 163 | c----------------------------------------------------------------------- |
|---|
| 164 | c longitude(imd) --> rlonud(imdp1) |
|---|
| 165 | c----------------------------------------------------------------------- |
|---|
| 166 | |
|---|
| 167 | c Passage en coordonnees boites scalaires et en radian |
|---|
| 168 | do i=1,imd |
|---|
| 169 | rlonud(i)=(longitude(i)+.5)*pi/180. |
|---|
| 170 | enddo |
|---|
| 171 | |
|---|
| 172 | c Repetition de la valeur im+1 |
|---|
| 173 | rlonud(imdp1)=rlonud(1) + 2*pi |
|---|
| 174 | |
|---|
| 175 | c----------------------------------------------------------------------- |
|---|
| 176 | c latitude(jmdp1) --> rlonvd(jmd) |
|---|
| 177 | c----------------------------------------------------------------------- |
|---|
| 178 | |
|---|
| 179 | c Passage en coordonnees boites scalaires et en radian |
|---|
| 180 | do j=1,jmd |
|---|
| 181 | rlatvd(j)=(latitude(j)-.5)*pi/180. |
|---|
| 182 | enddo |
|---|
| 183 | |
|---|
| 184 | c======================================================================= |
|---|
| 185 | c lecture NetCDF de albedo, thermal, relief, zdtm (pour francois Lott) |
|---|
| 186 | c======================================================================= |
|---|
| 187 | |
|---|
| 188 | string(0) = 'z0' |
|---|
| 189 | string(1) = 'albedo' |
|---|
| 190 | string(2) = 'thermal' |
|---|
| 191 | if (relief.ne.'pla') then |
|---|
| 192 | write(*,*) ' MOLA topography' |
|---|
| 193 | relief = 'MOL' |
|---|
| 194 | string(3) = 'z'//relief |
|---|
| 195 | else |
|---|
| 196 | string(3) = 'zMOL' ! pour qu''il lise qqchose sur le fichier |
|---|
| 197 | ! remise a 0 derriere |
|---|
| 198 | endif |
|---|
| 199 | string(4) = 'zMOL' ! lecture pour calcul topog. sous-maille |
|---|
| 200 | |
|---|
| 201 | |
|---|
| 202 | DO k=0,4 |
|---|
| 203 | write(*,*) 'string',k,string(k) |
|---|
| 204 | |
|---|
| 205 | c----------------------------------------------------------------------- |
|---|
| 206 | c initialisation |
|---|
| 207 | c----------------------------------------------------------------------- |
|---|
| 208 | call initial0(iimp1*jjp1,pfield) |
|---|
| 209 | call initial0(imd*jmdp1,zdata) |
|---|
| 210 | call initial0(imdp1*jmdp1,zdataS) |
|---|
| 211 | |
|---|
| 212 | c----------------------------------------------------------------------- |
|---|
| 213 | c Lecture NetCDF |
|---|
| 214 | c----------------------------------------------------------------------- |
|---|
| 215 | |
|---|
| 216 | ierr = NF_INQ_VARID (unit, string(k), nvarid) |
|---|
| 217 | if (ierr.ne.nf_noerr) then |
|---|
| 218 | write(*,*) 'datareadnc error, cannot find ',trim(string(k)) |
|---|
| 219 | write(*,*) ' in file ',trim(datadir),'/surface.nc' |
|---|
| 220 | stop |
|---|
| 221 | endif |
|---|
| 222 | #ifdef NC_DOUBLE |
|---|
| 223 | ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata) |
|---|
| 224 | #else |
|---|
| 225 | ierr = NF_GET_VAR_REAL(unit, nvarid, zdata) |
|---|
| 226 | #endif |
|---|
| 227 | if (ierr.ne.nf_noerr) then |
|---|
| 228 | write(*,*) 'datareadnc: error failed loading ',trim(string(k)) |
|---|
| 229 | stop |
|---|
| 230 | endif |
|---|
| 231 | |
|---|
| 232 | c----------------------------------------------------------------------- |
|---|
| 233 | c Cas particulier "Francois Lott" ( k=4 ) (relief sous-maille) |
|---|
| 234 | c----------------------------------------------------------------------- |
|---|
| 235 | if (k.eq.4) then |
|---|
| 236 | |
|---|
| 237 | zdata(:)=1000.*zdata(:) |
|---|
| 238 | longitude(:)=(pi/180.)*longitude(:) |
|---|
| 239 | latitude(:)=(pi/180.)*latitude(:) |
|---|
| 240 | |
|---|
| 241 | call grid_noro1(360, 180, longitude, latitude, zdata, |
|---|
| 242 | . iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe) |
|---|
| 243 | |
|---|
| 244 | !CW17 |
|---|
| 245 | call avg_horiz(imd,jmdp1,iim,jjm,longitude, |
|---|
| 246 | . latitude,rlonv,rlatu,zdata,pfield) |
|---|
| 247 | |
|---|
| 248 | do j=1,jjp1 !CW17 49, iimp1=65, the last column = first column |
|---|
| 249 | pfield(iimp1*j) = pfield(1+iimp1*(j-1)) |
|---|
| 250 | enddo |
|---|
| 251 | do i=1,iimp1*jjp1 |
|---|
| 252 | c if (pfield(i) .ne. -999999.) then |
|---|
| 253 | zavg(i) = pfield(i) |
|---|
| 254 | c else |
|---|
| 255 | c zavg(i)=-999999. |
|---|
| 256 | c endif |
|---|
| 257 | enddo |
|---|
| 258 | |
|---|
| 259 | endif |
|---|
| 260 | |
|---|
| 261 | c----------------------------------------------------------------------- |
|---|
| 262 | c Passage de zdata en grille (imdp1 jmdp1) |
|---|
| 263 | c----------------------------------------------------------------------- |
|---|
| 264 | do j=1,jmdp1 |
|---|
| 265 | do i=1,imd |
|---|
| 266 | zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1)) |
|---|
| 267 | enddo |
|---|
| 268 | zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1)) |
|---|
| 269 | enddo |
|---|
| 270 | |
|---|
| 271 | c----------------------------------------------------------------------- |
|---|
| 272 | c Interpolation |
|---|
| 273 | c----------------------------------------------------------------------- |
|---|
| 274 | call interp_horiz(zdataS,pfield,imd,jmd, |
|---|
| 275 | . iim, jjm,1,rlonud,rlatvd,rlonu,rlatv) |
|---|
| 276 | |
|---|
| 277 | c----------------------------------------------------------------------- |
|---|
| 278 | c Periodicite |
|---|
| 279 | c----------------------------------------------------------------------- |
|---|
| 280 | |
|---|
| 281 | do j=1,jjp1 |
|---|
| 282 | pfield(iimp1*j) = pfield(1+iimp1*(j-1)) |
|---|
| 283 | enddo |
|---|
| 284 | |
|---|
| 285 | c----------------------------------------------------------------------- |
|---|
| 286 | c Sauvegarde des champs |
|---|
| 287 | c----------------------------------------------------------------------- |
|---|
| 288 | |
|---|
| 289 | if (k.eq.0) then ! z0 |
|---|
| 290 | z0(1:iimp1*jjp1)=pfield(1:iimp1*jjp1)*.01 |
|---|
| 291 | ! multiplied by 0.01 to have z0 in m |
|---|
| 292 | elseif (k.eq.1) then ! albedo |
|---|
| 293 | do i=1,iimp1*jjp1 |
|---|
| 294 | alb(i) = pfield(i) |
|---|
| 295 | enddo |
|---|
| 296 | elseif (k.eq.2) then ! thermal |
|---|
| 297 | do i=1,iimp1*jjp1 |
|---|
| 298 | ith(i) = pfield(i) |
|---|
| 299 | enddo |
|---|
| 300 | elseif (k.eq.3) then ! relief |
|---|
| 301 | if (relief.eq.'pla') then |
|---|
| 302 | call initial0(iimp1*jjp1,phisinit) |
|---|
| 303 | else |
|---|
| 304 | do i=1,iimp1*jjp1 |
|---|
| 305 | phisinit(i) = pfield(i) |
|---|
| 306 | enddo |
|---|
| 307 | endif |
|---|
| 308 | endif |
|---|
| 309 | |
|---|
| 310 | ENDDO |
|---|
| 311 | |
|---|
| 312 | c----------------------------------------------------------------------- |
|---|
| 313 | c Traitement Phisinit |
|---|
| 314 | c----------------------------------------------------------------------- |
|---|
| 315 | |
|---|
| 316 | phisinit(1:iimp1*jjp1)=1000.*phisinit(1:iimp1*jjp1) |
|---|
| 317 | phisinit(:)=g*phisinit(:) |
|---|
| 318 | |
|---|
| 319 | c----------------------------------------------------------------------- |
|---|
| 320 | c FIN |
|---|
| 321 | c----------------------------------------------------------------------- |
|---|
| 322 | |
|---|
| 323 | c======================================================================= |
|---|
| 324 | c<<<<<<<<<<<<<<<<<<<<<<<add new dataset>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
|---|
| 325 | ! name of dataset |
|---|
| 326 | string(5) = 'hmons' !subgrid hmons |
|---|
| 327 | c the following data could be useful in future, but currently, |
|---|
| 328 | c we don't need to put them into restartfi.nc |
|---|
| 329 | string(6) = 'summit' !subgrid summit |
|---|
| 330 | string(7) = 'base' !base of summit (not subgrid) |
|---|
| 331 | c string(8) = 'bottom' !subgrid base |
|---|
| 332 | do k=5,7 |
|---|
| 333 | write(*,*) 'string',k,string(k) |
|---|
| 334 | c----------------------------------------------------------------------- |
|---|
| 335 | c Lecture NetCDF |
|---|
| 336 | c----------------------------------------------------------------------- |
|---|
| 337 | |
|---|
| 338 | ierr = NF_INQ_VARID (unit, string(k), nvarid) |
|---|
| 339 | if (ierr.ne.nf_noerr) then |
|---|
| 340 | write(*,*) 'datareadnc error, cannot find ',trim(string(k)) |
|---|
| 341 | write(*,*) ' in file ',trim(datadir),'/surface.nc' |
|---|
| 342 | stop |
|---|
| 343 | endif |
|---|
| 344 | |
|---|
| 345 | c----------------------------------------------------------------------- |
|---|
| 346 | c initialisation |
|---|
| 347 | c----------------------------------------------------------------------- |
|---|
| 348 | call initial0(iimp1*jjp1,pfield) |
|---|
| 349 | call initial0(imd*jmdp1,zdata) |
|---|
| 350 | call initial0(imdp1*jmdp1,zdataS) |
|---|
| 351 | |
|---|
| 352 | #ifdef NC_DOUBLE |
|---|
| 353 | ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata) |
|---|
| 354 | #else |
|---|
| 355 | ierr = NF_GET_VAR_REAL(unit, nvarid, zdata) |
|---|
| 356 | #endif |
|---|
| 357 | if (ierr.ne.nf_noerr) then |
|---|
| 358 | write(*,*) 'datareadnc: error failed loading ', |
|---|
| 359 | . trim(string(k)) |
|---|
| 360 | stop |
|---|
| 361 | endif |
|---|
| 362 | |
|---|
| 363 | c----------------------------------------------------------------------- |
|---|
| 364 | c Passage de zdata en grille (imdp1 jmdp1) |
|---|
| 365 | c----------------------------------------------------------------------- |
|---|
| 366 | do j=1,jmdp1 ! 180 |
|---|
| 367 | do i=1,imd ! 360 |
|---|
| 368 | !copy zdata to zdataS, line by line |
|---|
| 369 | ! i+ 361 *(j-1) i+ 360*(j-1) |
|---|
| 370 | zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1)) |
|---|
| 371 | enddo |
|---|
| 372 | ! the last column = the first column of zdata |
|---|
| 373 | zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1)) |
|---|
| 374 | enddo |
|---|
| 375 | |
|---|
| 376 | if (k .eq. 5 .or. k .eq. 6 .or. k .eq. 7) then |
|---|
| 377 | ! hmons, summit, keep the maximum of subgrids |
|---|
| 378 | call mvc_horiz(imd,jmdp1,iim,jjm,longitude, |
|---|
| 379 | . latitude,rlonv,rlatu,zdata,pfield) |
|---|
| 380 | endif |
|---|
| 381 | |
|---|
| 382 | |
|---|
| 383 | c----------------------------------------------------------------------- |
|---|
| 384 | c Periodicite |
|---|
| 385 | c----------------------------------------------------------------------- |
|---|
| 386 | |
|---|
| 387 | do j=1,jjp1 ! 49, iimp1=65, the last column = first column |
|---|
| 388 | pfield(iimp1*j) = pfield(1+iimp1*(j-1)) |
|---|
| 389 | enddo |
|---|
| 390 | |
|---|
| 391 | c----------------------------------------------------------------------- |
|---|
| 392 | c Sauvegarde des champs |
|---|
| 393 | c----------------------------------------------------------------------- |
|---|
| 394 | |
|---|
| 395 | if (k.eq.5) then ! hmons |
|---|
| 396 | do i=1,iimp1*jjp1 |
|---|
| 397 | if (pfield(i) .ne. -999999.) then |
|---|
| 398 | hmons(i) = pfield(i) |
|---|
| 399 | else |
|---|
| 400 | hmons(i)=0. |
|---|
| 401 | endif |
|---|
| 402 | enddo |
|---|
| 403 | endif |
|---|
| 404 | |
|---|
| 405 | if (k.eq.6) then ! summit |
|---|
| 406 | do i=1,iimp1*jjp1 |
|---|
| 407 | if (pfield(i) .ne. -999999.) then |
|---|
| 408 | summit(i) = pfield(i) |
|---|
| 409 | else |
|---|
| 410 | summit(i)=0. |
|---|
| 411 | endif |
|---|
| 412 | enddo |
|---|
| 413 | endif |
|---|
| 414 | |
|---|
| 415 | if (k.eq.7) then ! base |
|---|
| 416 | do i=1,iimp1*jjp1 |
|---|
| 417 | if (pfield(i) .ne. -999999.) then |
|---|
| 418 | base(i) = pfield(i) |
|---|
| 419 | else |
|---|
| 420 | base(i)=0. |
|---|
| 421 | endif |
|---|
| 422 | enddo |
|---|
| 423 | endif |
|---|
| 424 | |
|---|
| 425 | enddo |
|---|
| 426 | c<<<<<<<<<<<<<<<<<<<<<<<<<done add new dataset>>>>>>>>>>>>>>>>>>>>>>>>>> |
|---|
| 427 | c======================================================================= |
|---|
| 428 | END |
|---|
