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