source: trunk/LMDZ.MARS/libf/phymars/nirco2abs.F @ 2740

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

LMDZ_MARS RV : Open_MP;
Put all the "save" variables as "!$OMP THREADPRIVATE" in phymars.
The code can now be tested, see README for more info

File size: 9.1 KB
Line 
1      SUBROUTINE nirco2abs(ngrid,nlayer,pplay,dist_sol,nq,pq,
2     $     mu0,fract,declin,pdtnirco2)
3                                                   
4       use tracer_mod, only: igcm_co2, igcm_o
5       use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
6       USE comcstfi_h, ONLY: pi
7       USE time_phylmdz_mod, ONLY: daysec
8       IMPLICIT NONE
9c=======================================================================
10c   subject:
11c   --------
12c   Computing heating rate due to
13c   absorption by CO2 in the near-infrared
14c   This version includes NLTE effects
15c
16c   (Scheme to be described in Forget et al., JGR, 2003)
17c   (old Scheme described in Forget et al., JGR, 1999)
18c
19c   This version updated with a new functional fit,
20c   see NLTE correction-factor of Lopez-Valverde et al (1998)
21c   Stephen Lewis 2000
22c
23c   jul 2011 malv+fgg    New corrections for NLTE implemented
24c   08/2002 : correction for bug when running with diurnal=F
25c
26c   author:  Frederic Hourdin 1996
27c   ------
28c            Francois Forget 1999
29c
30c   input:
31c   -----
32c   ngrid                 number of gridpoint of horizontal grid
33c   nlayer                Number of layer
34c   dist_sol              sun-Mars distance (AU)
35c   mu0(ngrid)         
36c   fract(ngrid)          day fraction of the time interval
37c   declin                latitude of subslar point
38c
39c   output:
40c   -------
41c
42c   pdtnirco2(ngrid,nlayer)      Heating rate (K/s)
43c
44c
45c=======================================================================
46c
47c    0.  Declarations :
48c    ------------------
49c
50      include "callkeys.h"
51      include "nirdata.h"
52
53c-----------------------------------------------------------------------
54c    Input/Output
55c    ------------
56      integer,intent(in) :: ngrid ! number of (horizontal) grid points
57      integer,intent(in) :: nlayer ! number of atmospheric layers
58      real,intent(in) :: pplay(ngrid,nlayer) ! Pressure
59      real,intent(in) :: dist_sol ! Sun-Mars distance (in AU)
60      integer,intent(in) :: nq ! number of tracers
61      real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers
62      real,intent(in) :: mu0(ngrid) ! solar angle
63      real,intent(in) :: fract(ngrid) ! day fraction of the time interval
64      real,intent(in) :: declin ! latitude of sub-solar point
65     
66      real,intent(out) :: pdtnirco2(ngrid,nlayer) ! heating rate (K/s)
67c
68c    Local variables :
69c    -----------------
70      INTEGER l,ig, n, nstep,i
71      REAL co2heat0, zmu(ngrid)
72
73c     special diurnal=F
74      real mu0_int(ngrid),fract_int(ngrid),zday_int
75      real ztim1,ztim2,ztim3,step
76
77c
78c   local saved variables
79c   ---------------------
80      logical,save :: firstcall=.true.
81      integer,save :: ico2=0 ! index of "co2" tracer
82      integer,save :: io=0 ! index of "o" tracer
83     
84!$OMP THREADPRIVATE(firstcall,ico2,io)
85     
86c     p0noonlte is a pressure below which non LTE effects are significant.
87c     REAL p0nonlte
88c     DATA p0nonlte/7.5e-3/
89c     SAVE p0nonlte
90
91c     parameters for CO2 heating fit
92      real,parameter :: n_a=1.1956475
93      real,parameter :: n_p0=0.0015888279
94      real,parameter :: n_b=1.9628251
95
96c     Variables added to implement NLTE correction factor (feb 2011)
97      real    pyy(nlayer)
98      real    cor1(nlayer),oldoco2(nlayer),alfa2(nlayer)
99      real    p2011,cociente1,merge
100      real    cor0,oco2gcm
101
102c----------------------------------------------------------------------
103
104c     Initialisation
105c     --------------
106      ! AS: OK firstcall absolute
107      if (firstcall) then
108        if (nircorr.eq.1) then
109          ! we will need co2 and o tracers
110          ico2=igcm_co2
111          if (ico2==0) then
112            write(*,*) "nirco2abs error: I need a CO2 tracer"
113            write(*,*) "     when running with nircorr==1"
114            call abort_physic("nirco2abs","need a CO2 tracer",1)
115          endif
116          io=igcm_o
117          if (io==0) then
118            write(*,*) "nirco2abs error: I need an O tracer"
119            write(*,*) "     when running with nircorr==1"
120            call abort_physic("nirco2abs","need an O tracer",1)
121          endif
122        endif
123        firstcall=.false.
124      endif
125
126
127c     co2heat is the heating by CO2 at 700Pa for a zero zenithal angle.
128      co2heat0=n_a*(1.52/dist_sol)**2/daysec
129
130c     Simple calcul for a given sun incident angle (if diurnal=T)
131c     --------------------------------------------
132
133      IF (diurnal) THEN 
134         do ig=1,ngrid
135            zmu(ig)=sqrt(1224.*mu0(ig)*mu0(ig)+1.)/35.
136
137            if(nircorr.eq.1) then
138               do l=1,nlayer
139                  pyy(l)=pplay(ig,l)
140               enddo
141
142               call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
143               call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
144               call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
145            endif
146
147            do l=1,nlayer
148!           Calculations for the O/CO2 correction
149               if(nircorr.eq.1) then
150                  cor0=1./(1.+n_p0/pplay(ig,l))**n_b
151                  if(pq(ig,l,ico2).gt.1.e-6) then                   
152                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
153                     ! handle the rare cases when pq(ig,l,io)<0
154                     if (pq(ig,l,io).lt.0) then
155                       write(*,*) "nirco2abs: warning ig=",ig," l=",l,
156     &                            " pq(ig,l,io)=",pq(ig,l,io)
157                       oco2gcm=1.e6
158                     endif
159                  else
160                     oco2gcm=1.e6
161                  endif
162                  cociente1=oco2gcm/oldoco2(l)
163                  merge=alog10(cociente1)*alfa2(l)+alog10(cor0)*
164     $                 (1.-alfa2(l))
165                  merge=10**merge
166                  p2011=sqrt(merge)*cor0
167               else if (nircorr.eq.0) then
168                  p2011=1.
169                  cor1(l)=1.
170               endif
171
172               if(fract(ig).gt.0.) pdtnirco2(ig,l)=
173     &             co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l))
174     &             /(1.+n_p0/pplay(ig,l))**n_b
175!           Corrections from tabulation
176     $              * cor1(l) * p2011
177c           OLD SCHEME (forget et al. 1999)
178c    s           co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l))
179c    s          / (1+p0nonlte/pplay(ig,l))
180           enddo
181         enddo
182         
183
184c     Averaging over diurnal cycle (if diurnal=F)
185c     -------------------------------------------
186c     NIR CO2 abs is slightly non linear. To remove the diurnal
187c     cycle, it is better to average the heating rate over 1 day rather
188c     than using the mean mu0 computed by mucorr in physiq.F (FF, 1998)
189
190      ELSE      ! if (.not.diurnal) then
191
192         nstep = 20   ! number of integration step /sol
193         mu0_int(1:ngrid) = 0.
194         ztim1 = 0.
195         do n=1,nstep
196            zday_int = (n-1)/float(nstep)
197            ztim2=COS(declin)*COS(2.*pi*(zday_int-.5))
198            ztim3=-COS(declin)*SIN(2.*pi*(zday_int-.5))
199            CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
200     s             ztim1,ztim2,ztim3,
201     s             mu0_int,fract_int)
202            do ig=1,ngrid               
203               zmu(ig)=sqrt(1224.*mu0_int(ig)*mu0_int(ig)+1.)/35.
204
205               if(nircorr.eq.1) then
206                  do l=1,nlayer
207                     pyy(l)=pplay(ig,l)
208                  enddo
209               call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
210               call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
211               call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
212               endif
213
214               do l=1,nlayer
215                  if(nircorr.eq.1) then
216                     cor0=1./(1.+n_p0/pplay(ig,l))**n_b
217                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
218                     cociente1=oco2gcm/oldoco2(l)
219                     merge=alog10(cociente1)*alfa2(l)+alog10(cor0)*
220     $                    (1.-alfa2(l))
221                     merge=10**merge
222                     p2011=sqrt(merge)*cor0
223                  else if (nircorr.eq.0) then
224                     p2011=1.
225                     cor1(l)=1.
226                  endif
227
228                  if(fract_int(ig).gt.0.) pdtnirco2(ig,l)=
229     &                 pdtnirco2(ig,l) + (1/float(nstep))*
230     &                 co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l))
231     &                 /(1.+n_p0/pplay(ig,l))**n_b
232!                      Corrections from tabulation
233     $                 * cor1(l) * p2011
234               enddo
235            enddo
236         end do
237      END IF
238
239      end
240
241
242     
243      subroutine interpnir(escout,p,nlayer,escin,pin,nl)
244C
245C subroutine to perform linear interpolation in pressure from 1D profile
246C escin(nl) sampled on pressure grid pin(nl) to profile
247C escout(nlayer) on pressure grid p(nlayer).
248C
249      real escout(nlayer),p(nlayer)
250      real escin(nl),pin(nl),wm,wp
251      integer nl,nlayer,n1,n,nm,np
252      do n1=1,nlayer
253         if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then
254            escout(n1) = 0.0
255         else
256            do n = 1,nl-1
257               if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
258                  nm=n
259                  np=n+1
260                  wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
261                  wp=1.0 - wm
262               endif
263            enddo
264            escout(n1) = escin(nm)*wm + escin(np)*wp
265         endif
266      enddo
267      end
Note: See TracBrowser for help on using the repository browser.