source: trunk/LMDZ.MARS/libf/phymars/read_dust_scenario.F90 @ 2616

Last change on this file since 2616 was 2586, checked in by romain.vande, 3 years ago

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