| 1 | subroutine read_dust_scenario(ngrid,nlayer,zday,pplev,tau_pref_scenario) |
|---|
| 2 | |
|---|
| 3 | ! Reading of the dust scenario file |
|---|
| 4 | |
|---|
| 5 | use netcdf |
|---|
| 6 | use geometry_mod, only: latitude, longitude ! in radians |
|---|
| 7 | use datafile_mod, only: datadir |
|---|
| 8 | use dust_param_mod, only: odpref, dustscaling_mode |
|---|
| 9 | use planete_h, only: year_day |
|---|
| 10 | implicit none |
|---|
| 11 | |
|---|
| 12 | include "callkeys.h" |
|---|
| 13 | |
|---|
| 14 | integer, intent(in) :: ngrid,nlayer |
|---|
| 15 | real, intent(in) :: zday ! date in martian days |
|---|
| 16 | real, dimension(ngrid,nlayer+1), intent(in) :: pplev |
|---|
| 17 | real, dimension(ngrid), intent(out) :: tau_pref_scenario ! visible dust column |
|---|
| 18 | ! opacity at odpref from scenario |
|---|
| 19 | |
|---|
| 20 | ! Local variables |
|---|
| 21 | |
|---|
| 22 | real :: realday |
|---|
| 23 | integer nid,nvarid,ierr |
|---|
| 24 | integer ltloop,lsloop,iloop,jloop,varloop,ig |
|---|
| 25 | real, dimension(2) :: taubuf |
|---|
| 26 | real tau1(4),tau |
|---|
| 27 | real alt(4) |
|---|
| 28 | integer latp(4),lonp(4) |
|---|
| 29 | integer yinf,ysup,xinf,xsup,tinf,tsup |
|---|
| 30 | real latinf,latsup,loninf,lonsup |
|---|
| 31 | real latintmp,lonintmp,latdeg,londeg |
|---|
| 32 | real colat,dlat,dlon,colattmp |
|---|
| 33 | logical, save :: firstcall=.true. |
|---|
| 34 | logical :: timeflag |
|---|
| 35 | real,save :: radeg,pi |
|---|
| 36 | integer :: timedim,londim,latdim |
|---|
| 37 | real, dimension(:), allocatable, save :: lat,lon,time |
|---|
| 38 | real, dimension(:,:,:), allocatable, save :: tautes |
|---|
| 39 | integer, save :: timelen,lonlen,latlen |
|---|
| 40 | character(len=33),save :: filename |
|---|
| 41 | |
|---|
| 42 | realday=mod(zday,669.) |
|---|
| 43 | |
|---|
| 44 | ! AS: firstcall OK absolute |
|---|
| 45 | if (firstcall) then |
|---|
| 46 | firstcall=.false. |
|---|
| 47 | |
|---|
| 48 | pi=acos(-1.) |
|---|
| 49 | radeg=180/pi |
|---|
| 50 | |
|---|
| 51 | ! assimilated dust file: (NB: iaervar is a common in "callkeys.h") |
|---|
| 52 | ! iaervar=4 means read dust_tes.nc file |
|---|
| 53 | ! iaervar=6 means read dust_cold.nc file |
|---|
| 54 | ! iaervar=7 means read dust_warm.nc file |
|---|
| 55 | ! iaervar=8 means read dust_clim.nc file |
|---|
| 56 | ! iaervar=24 means read dust_MY24.nc file |
|---|
| 57 | ! iaervar=25 means read dust_MY25.nc file |
|---|
| 58 | ! iaervar=26 means read dust_MY26.nc file, etc. |
|---|
| 59 | if (iaervar.eq.4) then |
|---|
| 60 | ! NB: 4: old TES assimilated MY24 dust scenarios (at 700Pa ref pressure!) |
|---|
| 61 | filename="dust_tes.nc" |
|---|
| 62 | else if (iaervar.eq.6) then |
|---|
| 63 | filename="dust_cold.nc" |
|---|
| 64 | else if (iaervar.eq.7) then |
|---|
| 65 | filename="dust_warm.nc" |
|---|
| 66 | else if (iaervar.eq.8) then |
|---|
| 67 | filename="dust_clim.nc" |
|---|
| 68 | else if ((iaervar.ge.24).and.(iaervar.le.35)) then |
|---|
| 69 | write(filename,'(a7,i2,a3)')"dust_MY",iaervar,".nc" |
|---|
| 70 | ! 124,125,126: old TES assimilated dust scenarios (at 700Pa ref pressure!) |
|---|
| 71 | else if (iaervar.eq.124) then |
|---|
| 72 | filename="dust_tes_MY24.nc" |
|---|
| 73 | else if (iaervar.eq.125) then |
|---|
| 74 | filename="dust_tes_MY25.nc" |
|---|
| 75 | else if (iaervar.eq.126) then |
|---|
| 76 | filename="dust_tes_MY26.nc" |
|---|
| 77 | endif |
|---|
| 78 | |
|---|
| 79 | ierr=nf90_open(trim(datadir)//"/"//trim(filename),NF90_NOWRITE,nid) |
|---|
| 80 | IF (ierr.NE.nf90_noerr) THEN |
|---|
| 81 | write(*,*)'Problem opening ',trim(filename),' (in phymars/read_dust_scenario.F90)' |
|---|
| 82 | write(*,*)'It should be in :',trim(datadir),'/' |
|---|
| 83 | write(*,*)'1) You can change this directory address in callfis.def with' |
|---|
| 84 | write(*,*)' datadir=/path/to/datafiles' |
|---|
| 85 | write(*,*)'2) If necessary, dust*.nc (and other datafiles)' |
|---|
| 86 | write(*,*)' can be obtained online on:' |
|---|
| 87 | write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir' |
|---|
| 88 | CALL abort_physic("read_dust_scenario","Cannot find file",1) |
|---|
| 89 | ENDIF |
|---|
| 90 | |
|---|
| 91 | ierr=nf90_inq_dimid(nid,"Time",timedim) |
|---|
| 92 | if (ierr.ne.nf90_noerr) then |
|---|
| 93 | ! 'Time' dimension not found, look for 'time' |
|---|
| 94 | ierr=nf90_inq_dimid(nid,"time",timedim) |
|---|
| 95 | if (ierr.ne.nf90_noerr) then |
|---|
| 96 | write(*,*)"Error: read_dust_scenario <time> not found" |
|---|
| 97 | endif |
|---|
| 98 | endif |
|---|
| 99 | ierr=nf90_inquire_dimension(nid,timedim,len=timelen) |
|---|
| 100 | ierr=nf90_inq_dimid(nid,"latitude",latdim) |
|---|
| 101 | ierr=nf90_inquire_dimension(nid,latdim,len=latlen) |
|---|
| 102 | ierr=nf90_inq_dimid(nid,"longitude",londim) |
|---|
| 103 | ierr=nf90_inquire_dimension(nid,londim,len=lonlen) |
|---|
| 104 | |
|---|
| 105 | |
|---|
| 106 | if (.not.allocated(tautes)) allocate(tautes(lonlen,latlen,timelen)) |
|---|
| 107 | if (.not.allocated(lat)) allocate(lat(latlen), lon(lonlen), time(timelen)) |
|---|
| 108 | |
|---|
| 109 | ! "dustop" if loading visible extinction opacity |
|---|
| 110 | ! "cdod" if loading IR absorption opacity |
|---|
| 111 | ierr=nf90_inq_varid(nid,"dustop",nvarid) |
|---|
| 112 | if (ierr.eq.nf90_noerr) then |
|---|
| 113 | ierr=nf90_get_var(nid,nvarid,tautes) |
|---|
| 114 | IF (ierr .NE. nf90_noerr) THEN |
|---|
| 115 | PRINT*, "Error: read_dust_scenario <dustop> not found" |
|---|
| 116 | write(*,*)trim(nf90_strerror(ierr)) |
|---|
| 117 | call abort_physic("read_dust_scenario","dustop not found",1) |
|---|
| 118 | ENDIF |
|---|
| 119 | else |
|---|
| 120 | ! did not find "dustop" , look for "cdod" |
|---|
| 121 | ierr=nf90_inq_varid(nid,"cdod",nvarid) |
|---|
| 122 | ierr=nf90_get_var(nid,nvarid,tautes) |
|---|
| 123 | IF (ierr .NE. nf90_noerr) THEN |
|---|
| 124 | PRINT*, "Error: read_dust_scenario <cdod> not found" |
|---|
| 125 | write(*,*)trim(nf90_strerror(ierr)) |
|---|
| 126 | call abort_physic("read_dust_scenario","cdod not found",1) |
|---|
| 127 | ENDIF |
|---|
| 128 | ! and multiply by 2*1.3=2.6 to convert from IR absorption |
|---|
| 129 | ! to visible extinction opacity |
|---|
| 130 | tautes(:,:,:)=2.6*tautes(:,:,:) |
|---|
| 131 | endif |
|---|
| 132 | |
|---|
| 133 | ierr=nf90_inq_varid(nid,"Time",nvarid) |
|---|
| 134 | ierr=nf90_get_var(nid,nvarid,time) |
|---|
| 135 | IF (ierr .NE. nf90_noerr) THEN |
|---|
| 136 | PRINT*, "Error: read_dust_scenario <Time> not found" |
|---|
| 137 | write(*,*)trim(nf90_strerror(ierr)) |
|---|
| 138 | call abort_physic("read_dust_scenario","Time not found",1) |
|---|
| 139 | ENDIF |
|---|
| 140 | |
|---|
| 141 | ! Adapt time axis values if using dustinjection scheme and/or |
|---|
| 142 | ! dustscaling_mode == 2 |
|---|
| 143 | IF ((dustinjection>0).or.(dustscaling_mode==2)) THEN |
|---|
| 144 | ! Sanity check: we expect to have one map per sol |
|---|
| 145 | if (timelen/=year_day) then |
|---|
| 146 | write(*,*)"read_dust_scenario error: timelen=",timelen |
|---|
| 147 | write(*,*)" whereas it was exepceted to be year_day=",year_day |
|---|
| 148 | call abort_physic("read_dust_scenario","timelen/=year_day",1) |
|---|
| 149 | endif |
|---|
| 150 | ! fill time() with intergers from 0 to timelen-1 |
|---|
| 151 | write(*,*)"read_dust_scenario: adapting time axis to integer values" |
|---|
| 152 | do iloop=1,timelen |
|---|
| 153 | time(iloop)=iloop-1 |
|---|
| 154 | enddo |
|---|
| 155 | ENDIF ! of IF ((dustinjection>0).or.(dustscaling_mode==2)) |
|---|
| 156 | |
|---|
| 157 | ierr=nf90_inq_varid(nid,"latitude",nvarid) |
|---|
| 158 | ierr=nf90_get_var(nid,nvarid,lat) |
|---|
| 159 | IF (ierr .NE. nf90_noerr) THEN |
|---|
| 160 | PRINT*, "Error: read_dust_scenario <latitude> not found" |
|---|
| 161 | write(*,*)trim(nf90_strerror(ierr)) |
|---|
| 162 | call abort_physic("read_dust_scenario","latitude not found",1) |
|---|
| 163 | ENDIF |
|---|
| 164 | |
|---|
| 165 | ierr=nf90_inq_varid(nid,"longitude",nvarid) |
|---|
| 166 | ierr=nf90_get_var(nid,nvarid,lon) |
|---|
| 167 | IF (ierr .NE. nf90_noerr) THEN |
|---|
| 168 | PRINT*, "Error: read_dust_scenario <longitude> not found" |
|---|
| 169 | write(*,*)trim(nf90_strerror(ierr)) |
|---|
| 170 | call abort_physic("read_dust_scenario","longitude not found",1) |
|---|
| 171 | ENDIF |
|---|
| 172 | |
|---|
| 173 | ierr=nf90_close(nid) |
|---|
| 174 | |
|---|
| 175 | endif ! of if (firstcall) |
|---|
| 176 | |
|---|
| 177 | do ig=1,ngrid |
|---|
| 178 | |
|---|
| 179 | ! Find the four nearest points, arranged as follows: |
|---|
| 180 | ! 1 2 |
|---|
| 181 | ! 3 4 |
|---|
| 182 | |
|---|
| 183 | latdeg=latitude(ig)*radeg ! latitude, in degrees |
|---|
| 184 | londeg=longitude(ig)*radeg ! longitude, in degrees east |
|---|
| 185 | colat=90-latdeg ! colatitude, in degrees |
|---|
| 186 | |
|---|
| 187 | ! Ehouarn: rounding effects and/or specific compiler issues |
|---|
| 188 | ! sometime cause londeg to be sligthly below -180 ... |
|---|
| 189 | if (londeg.lt.-180) then |
|---|
| 190 | ! check if it is by a large amount |
|---|
| 191 | if ((londeg+180).lt.-1.e-3) then |
|---|
| 192 | write(*,*) ' ig=',ig,' londeg=',londeg |
|---|
| 193 | call abort_physic("read_dust_scenario","longitude way out of range",1) |
|---|
| 194 | else |
|---|
| 195 | londeg=-180 |
|---|
| 196 | endif |
|---|
| 197 | endif |
|---|
| 198 | |
|---|
| 199 | ! Find encompassing latitudes |
|---|
| 200 | if (colat<(90-lat(1))) then ! between north pole and lat(1) |
|---|
| 201 | ysup=1 |
|---|
| 202 | yinf=1 |
|---|
| 203 | else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten) |
|---|
| 204 | ysup=latlen |
|---|
| 205 | yinf=latlen |
|---|
| 206 | else ! general case |
|---|
| 207 | do iloop=2,latlen |
|---|
| 208 | if(colat<(90-lat(iloop))) then |
|---|
| 209 | ysup=iloop-1 |
|---|
| 210 | yinf=iloop |
|---|
| 211 | exit |
|---|
| 212 | endif |
|---|
| 213 | enddo |
|---|
| 214 | endif |
|---|
| 215 | |
|---|
| 216 | latinf=lat(yinf) |
|---|
| 217 | latsup=lat(ysup) |
|---|
| 218 | |
|---|
| 219 | |
|---|
| 220 | ! Find encompassing longitudes |
|---|
| 221 | ! Note: in input file, lon(1)=-180. |
|---|
| 222 | if (londeg>lon(lonlen)) then |
|---|
| 223 | xsup=1 |
|---|
| 224 | xinf=lonlen |
|---|
| 225 | loninf=lon(xsup) |
|---|
| 226 | lonsup=180.0 ! since lon(1)=-180.0 |
|---|
| 227 | else |
|---|
| 228 | do iloop=2,lonlen |
|---|
| 229 | if(londeg<lon(iloop)) then |
|---|
| 230 | xsup=iloop |
|---|
| 231 | xinf=iloop-1 |
|---|
| 232 | exit |
|---|
| 233 | endif |
|---|
| 234 | enddo |
|---|
| 235 | loninf=lon(xinf) |
|---|
| 236 | lonsup=lon(xsup) |
|---|
| 237 | endif |
|---|
| 238 | |
|---|
| 239 | if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)& |
|---|
| 240 | .OR.(ysup.lt.1)) then |
|---|
| 241 | write (*,*) "read_dust_scenario: SYSTEM ERROR on x or y in ts_gcm" |
|---|
| 242 | write (*,*) "xinf: ",xinf |
|---|
| 243 | write (*,*) "xsup: ",xsup |
|---|
| 244 | write (*,*) "yinf: ",yinf |
|---|
| 245 | write (*,*) "ysup: ",ysup |
|---|
| 246 | call abort_physic("read_dust_scenario","invalid coordinates",1) |
|---|
| 247 | endif |
|---|
| 248 | |
|---|
| 249 | ! loninf=lon(xinf) |
|---|
| 250 | ! lonsup=lon(xsup) |
|---|
| 251 | ! latinf=lat(yinf) |
|---|
| 252 | ! latsup=lat(ysup) |
|---|
| 253 | |
|---|
| 254 | ! The four neighbouring points are arranged as follows: |
|---|
| 255 | ! 1 2 |
|---|
| 256 | ! 3 4 |
|---|
| 257 | |
|---|
| 258 | latp(1)=ysup |
|---|
| 259 | latp(2)=ysup |
|---|
| 260 | latp(3)=yinf |
|---|
| 261 | latp(4)=yinf |
|---|
| 262 | |
|---|
| 263 | lonp(1)=xinf |
|---|
| 264 | lonp(2)=xsup |
|---|
| 265 | lonp(3)=xinf |
|---|
| 266 | lonp(4)=xsup |
|---|
| 267 | |
|---|
| 268 | ! Linear interpolation on time, for all four neighbouring points |
|---|
| 269 | |
|---|
| 270 | if ((realday<time(1)).or.(realday>time(timelen))) then |
|---|
| 271 | tinf=timelen |
|---|
| 272 | tsup=1 |
|---|
| 273 | timeflag=.true. |
|---|
| 274 | else |
|---|
| 275 | timeflag=.false. |
|---|
| 276 | do iloop=2,timelen |
|---|
| 277 | if (realday<=time(iloop)) then |
|---|
| 278 | tinf=iloop-1 |
|---|
| 279 | tsup=iloop |
|---|
| 280 | exit |
|---|
| 281 | endif |
|---|
| 282 | enddo |
|---|
| 283 | endif |
|---|
| 284 | |
|---|
| 285 | ! Bilinear interpolation on the four nearest points |
|---|
| 286 | |
|---|
| 287 | if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then |
|---|
| 288 | dlat=0 |
|---|
| 289 | else |
|---|
| 290 | dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup)) |
|---|
| 291 | endif |
|---|
| 292 | |
|---|
| 293 | if (lonsup==loninf) then |
|---|
| 294 | dlon=0 |
|---|
| 295 | else |
|---|
| 296 | dlon=(londeg-loninf)/(lonsup-loninf) |
|---|
| 297 | endif |
|---|
| 298 | |
|---|
| 299 | do iloop=1,4 |
|---|
| 300 | taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf) |
|---|
| 301 | taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup) |
|---|
| 302 | if (timeflag) then |
|---|
| 303 | if (realday>time(timelen)) then |
|---|
| 304 | tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf))) |
|---|
| 305 | else |
|---|
| 306 | tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf))) |
|---|
| 307 | endif |
|---|
| 308 | else |
|---|
| 309 | tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf)) |
|---|
| 310 | endif |
|---|
| 311 | if (tau1(iloop)<0) then |
|---|
| 312 | write (*,*) "read_dust_scenario: SYSTEM ERROR on tau" |
|---|
| 313 | write (*,*) "utime ",realday |
|---|
| 314 | write (*,*) "time(tinf) ",time(tinf) |
|---|
| 315 | write (*,*) "time(tsup) ",time(tsup) |
|---|
| 316 | write (*,*) "tau1 ",taubuf(1) |
|---|
| 317 | write (*,*) "tau2 ",taubuf(2) |
|---|
| 318 | write (*,*) "tau ",tau1(iloop) |
|---|
| 319 | call abort_physic("read_dust_scenario","negative tau!",1) |
|---|
| 320 | endif |
|---|
| 321 | enddo |
|---|
| 322 | |
|---|
| 323 | if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then |
|---|
| 324 | write (*,*) "read_dust_scenario: SYSTEM ERROR on dlat or dlon in ts_gcm" |
|---|
| 325 | write (*,*) "dlat: ",dlat |
|---|
| 326 | write (*,*) "lat: ",latdeg |
|---|
| 327 | write (*,*) "dlon: ",dlon |
|---|
| 328 | write (*,*) "lon: ",londeg |
|---|
| 329 | call abort_physic("read_dust_scenario","negative dlat or dlon!",1) |
|---|
| 330 | endif |
|---|
| 331 | |
|---|
| 332 | tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3) |
|---|
| 333 | |
|---|
| 334 | tau_pref_scenario(ig)=tau |
|---|
| 335 | ! |
|---|
| 336 | enddo ! of ig=1,ngrid |
|---|
| 337 | |
|---|
| 338 | if (filename(1:8)=="dust_tes") then |
|---|
| 339 | ! when using old TES files, correction for: |
|---|
| 340 | ! - the reference pressure was 700Pa (unlike 610Pa now) |
|---|
| 341 | ! - the 1.3 conversion factor from IR absorbtion opacity to |
|---|
| 342 | ! IR extinction opacity |
|---|
| 343 | tau_pref_scenario(:)=tau_pref_scenario(:)*1.3*(odpref/700.) |
|---|
| 344 | endif |
|---|
| 345 | |
|---|
| 346 | if (swrtype.eq.1) then ! Fouquart (NB: swrtype is set in callkeys.h) |
|---|
| 347 | ! when using old radiative transfer (like in MCD 4.x) |
|---|
| 348 | ! needed to decrease opacity (*0.825) to compensate overestimation of |
|---|
| 349 | ! heating rates |
|---|
| 350 | tau_pref_scenario(:)=tau_pref_scenario(:)*0.825/1.3 |
|---|
| 351 | endif |
|---|
| 352 | |
|---|
| 353 | end |
|---|
