Changeset 4198
- Timestamp:
- Dec 20, 2019, 12:00:02 PM (5 years ago)
- Location:
- dynamico_lmdz/simple_physics/phyparam
- Files:
-
- 2 deleted
- 2 edited
-
param/gather.F (deleted)
-
param/phyparam.F (modified) (1 diff)
-
param/sw.F (deleted)
-
physics/radiative.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
dynamico_lmdz/simple_physics/phyparam/param/phyparam.F
r4197 r4198 14 14 USE astronomy 15 15 USE vdif_mod, ONLY : vdif 16 USE radiative, ONLY : mucorr 16 USE radiative, ONLY : mucorr, sw 17 17 c 18 18 IMPLICIT NONE -
dynamico_lmdz/simple_physics/phyparam/physics/radiative.F90
r4197 r4198 2 2 IMPLICIT NONE 3 3 SAVE 4 5 PRIVATE 6 7 PUBLIC :: mucorr, sw 4 8 5 9 CONTAINS 6 7 SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad) 8 10 11 PURE SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad) 12 13 !======================================================================= 14 ! 15 ! Calcul of equivalent solar angle and and fraction of day whithout 16 ! diurnal cycle. 17 ! 18 ! parmeters : 19 ! ----------- 20 ! 21 ! Input : 22 ! ------- 23 ! npts number of points 24 ! pdeclin solar declinaison 25 ! plat(npts) latitude 26 ! phaut hauteur typique de l'atmosphere 27 ! prad rayon planetaire 28 ! 29 ! Output : 30 ! -------- 31 ! pmu(npts) equivalent cosinus of the solar angle 32 ! pfract(npts) fractionnal day 33 ! 34 !======================================================================= 35 36 !----------------------------------------------------------------------- 37 ! 38 ! 0. Declarations : 39 ! ----------------- 40 41 ! Arguments : 42 ! ----------- 43 INTEGER, INTENT(IN) :: npts 44 REAL, INTENT(IN) :: phaut, prad, pdeclin, plat(npts) 45 REAL, INTENT(OUT) :: pmu(npts), pfract(npts) 46 ! 47 ! Local variables : 48 ! ----------------- 49 INTEGER j 50 REAL z,cz,sz,tz,phi,cphi,sphi,tphi 51 REAL ap,a,t,b 52 REAL alph 53 54 REAL, PARAMETER :: pi=2.*asin(1.) 55 56 !----------------------------------------------------------------------- 57 58 ! print*,'npts,pdeclin',npts,pdeclin 59 z = pdeclin 60 cz = cos (z) 61 sz = sin (z) 62 ! print*,'cz,sz',cz,sz 63 64 DO j = 1, npts 65 66 phi = plat(j) 67 cphi = cos(phi) 68 if (cphi.le.1.e-9) cphi=1.e-9 69 sphi = sin(phi) 70 tphi = sphi / cphi 71 b = cphi * cz 72 t = -tphi * sz / cz 73 a = 1.0 - t*t 74 ap = a 75 76 IF(t.eq.0.) then 77 t=0.5*pi 78 ELSE 79 IF (a.lt.0.) a = 0. 80 t = sqrt(a) / t 81 IF (t.lt.0.) then 82 t = -atan (-t) + pi 83 ELSE 84 t = atan(t) 85 ENDIF 86 ENDIF 87 88 pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi 89 pfract(j) = t / pi 90 IF (ap .lt.0.) then 91 pmu(j) = sphi * sz 92 pfract(j) = 1.0 93 ENDIF 94 IF (pmu(j).le.0.0) pmu(j) = 0.0 95 pmu(j) = pmu(j) / pfract(j) 96 IF (pmu(j).eq.0.) pfract(j) = 0. 97 98 END DO 99 100 !----------------------------------------------------------------------- 101 ! correction de rotondite: 102 ! ------------------------ 103 104 alph=phaut/prad 105 DO j=1,npts 106 ! !!!!!! 107 pmu(j)=sqrt(1224.*pmu(j)*pmu(j)+1.)/35. 108 ! $ (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j)) 109 END DO 110 111 END SUBROUTINE mucorr 112 113 PURE SUBROUTINE monGATHER(n,a,b,index) 114 INTEGER, INTENT(IN) :: n,index(n) 115 REAL, INTENT(IN) :: b(n) 116 REAL, INTENT(OUT) :: a(n) 117 INTEGER :: i 118 119 DO i=1,n 120 a(i)=b(index(i)) 121 END DO 122 END SUBROUTINE monGATHER 123 124 SUBROUTINE sw(ngrid,nlayer,ldiurn, & 125 coefvis,albedo, & 126 plevel,ps_rad,pmu,pfract,psolarf0, & 127 fsrfvis,dtsw, & 128 lwrite) 129 USE phys_const 9 130 !======================================================================= 10 131 ! 11 ! Calcul of equivalent solar angle and and fraction of day whithout 12 ! diurnal cycle. 13 ! 14 ! parmeters : 15 ! ----------- 16 ! 17 ! Input : 18 ! ------- 19 ! npts number of points 20 ! pdeclin solar declinaison 21 ! plat(npts) latitude 22 ! phaut hauteur typique de l'atmosphere 23 ! prad rayon planetaire 24 ! 25 ! Output : 26 ! -------- 27 ! pmu(npts) equivalent cosinus of the solar angle 28 ! pfract(npts) fractionnal day 132 ! Rayonnement solaire en atmosphere non diffusante avec un 133 ! coefficient d'absoprption gris. 29 134 ! 30 135 !======================================================================= 31 32 !----------------------------------------------------------------------- 33 ! 34 ! 0. Declarations : 35 ! ----------------- 36 37 ! Arguments : 38 ! ----------- 39 INTEGER npts 40 REAL plat(npts),pmu(npts), pfract(npts) 41 REAL phaut,prad,pdeclin 42 ! 43 ! Local variables : 44 ! ----------------- 45 INTEGER j 46 REAL pi 47 REAL z,cz,sz,tz,phi,cphi,sphi,tphi 48 REAL ap,a,t,b 49 REAL alph 50 51 !----------------------------------------------------------------------- 52 53 print*,'npts,pdeclin' 54 print*,npts,pdeclin 55 pi = 4. * atan(1.0) 56 print*,'PI=',pi 57 pi=2.*asin(1.) 58 print*,'PI=B',pi 59 z = pdeclin 60 cz = cos (z) 61 sz = sin (z) 62 print*,'cz,sz',cz,sz 63 64 DO j = 1, npts 65 66 phi = plat(j) 67 cphi = cos(phi) 68 if (cphi.le.1.e-9) cphi=1.e-9 69 sphi = sin(phi) 70 tphi = sphi / cphi 71 b = cphi * cz 72 t = -tphi * sz / cz 73 a = 1.0 - t*t 74 ap = a 75 76 IF(t.eq.0.) then 77 t=0.5*pi 78 ELSE 79 IF (a.lt.0.) a = 0. 80 t = sqrt(a) / t 81 IF (t.lt.0.) then 82 t = -atan (-t) + pi 83 ELSE 84 t = atan(t) 136 ! 137 ! declarations: 138 ! ------------- 139 ! 140 ! 141 ! arguments: 142 ! ---------- 143 ! 144 INTEGER ngrid,nlayer 145 REAL albedo(ngrid),coefvis 146 REAL pmu(ngrid),pfract(ngrid) 147 REAL plevel(ngrid,nlayer+1),ps_rad 148 REAL psolarf0 149 REAL fsrfvis(ngrid),dtsw(ngrid,nlayer) 150 LOGICAL lwrite,ldiurn 151 ! 152 ! variables locales: 153 ! ------------------ 154 ! 155 156 REAL zalb(ngrid),zmu(ngrid),zfract(ngrid) 157 REAL zplev(ngrid,nlayer+1) 158 REAL zflux(ngrid),zdtsw(ngrid,nlayer) 159 160 INTEGER ig,l,nlevel,index(ngrid),ncount,igout 161 REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1) 162 REAL zfsrfref(ngrid) 163 REAL z1(ngrid) 164 REAL zu(ngrid,nlayer+1) 165 REAL tau0 166 167 !----------------------------------------------------------------------- 168 ! 1. initialisations: 169 ! ------------------- 170 171 172 nlevel=nlayer+1 173 174 !----------------------------------------------------------------------- 175 ! Definitions des tableaux locaux pour les points ensoleilles: 176 ! ------------------------------------------------------------ 177 178 IF (ldiurn) THEN 179 ncount=0 180 DO ig=1,ngrid 181 index(ig)=0 182 ENDDO 183 DO ig=1,ngrid 184 IF(pfract(ig).GT.1.e-6) THEN 185 ncount=ncount+1 186 index(ncount)=ig 85 187 ENDIF 86 ENDIF 87 88 pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi 89 pfract(j) = t / pi 90 IF (ap .lt.0.) then 91 pmu(j) = sphi * sz 92 pfract(j) = 1.0 93 ENDIF 94 IF (pmu(j).le.0.0) pmu(j) = 0.0 95 pmu(j) = pmu(j) / pfract(j) 96 IF (pmu(j).eq.0.) pfract(j) = 0. 97 98 END DO 99 100 !----------------------------------------------------------------------- 101 ! correction de rotondite: 188 ENDDO 189 CALL monGATHER(ncount,zfract,pfract,index) 190 CALL monGATHER(ncount,zmu,pmu,index) 191 CALL monGATHER(ncount,zalb,albedo,index) 192 DO l=1,nlevel 193 CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index) 194 ENDDO 195 ELSE 196 ncount=ngrid 197 zfract(:)=pfract(:) 198 zmu(:)=pmu(:) 199 zalb(:)=albedo(:) 200 zplev(:,:)=plevel(:,:) 201 ENDIF 202 203 !----------------------------------------------------------------------- 204 ! calcul des profondeurs optiques integres depuis p=0: 205 ! ---------------------------------------------------- 206 207 tau0=-.5*log(coefvis) 208 209 ! calcul de la partie homogene de l'opacite 210 tau0=tau0/ps_rad 211 DO l=1,nlayer+1 212 DO ig=1,ncount 213 zu(ig,l)=tau0*zplev(ig,l) 214 ENDDO 215 ENDDO 216 217 !----------------------------------------------------------------------- 218 ! 2. calcul de la transmission depuis le sommet de l'atmosphere: 219 ! ----------------------------------------------------------- 220 221 DO l=1,nlevel 222 DO ig=1,ncount 223 ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig)) 224 ENDDO 225 ENDDO 226 227 IF (lwrite) THEN 228 igout=ncount/2+1 229 PRINT* 230 PRINT*,'Diagnostique des transmission dans le spectre solaire' 231 PRINT*,'zfract, zmu, zalb' 232 PRINT*,zfract(igout), zmu(igout), zalb(igout) 233 PRINT*,'Pression, quantite d abs, transmission' 234 DO l=1,nlayer+1 235 PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l) 236 ENDDO 237 ENDIF 238 239 !----------------------------------------------------------------------- 240 ! 3. taux de chauffage, ray. solaire direct: 241 ! ------------------------------------------ 242 243 DO l=1,nlayer 244 DO ig=1,ncount 245 zdtsw(ig,l)=g*psolarf0*zfract(ig)*zmu(ig)* & 246 (ztrdir(ig,l+1)-ztrdir(ig,l))/ & 247 (cpp*(zplev(ig,l)-zplev(ig,l+1))) 248 ENDDO 249 ENDDO 250 IF (lwrite) THEN 251 PRINT* 252 PRINT*,'Diagnostique des taux de chauffage solaires:' 253 PRINT*,' 1 taux de chauffage lie au ray. solaire direct' 254 DO l=1,nlayer 255 PRINT*,zdtsw(igout,l) 256 ENDDO 257 ENDIF 258 259 260 !----------------------------------------------------------------------- 261 ! 4. calcul du flux solaire arrivant sur le sol: 262 ! ---------------------------------------------- 263 264 DO ig=1,ncount 265 z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1) 266 zflux(ig)=(1.-zalb(ig))*z1(ig) 267 zfsrfref(ig)= zalb(ig)*z1(ig) 268 ENDDO 269 IF (lwrite) THEN 270 PRINT* 271 PRINT*,'Diagnostique des taux de chauffage solaires:' 272 PRINT*,' 2 flux solaire net incident sur le sol' 273 PRINT*,zflux(igout) 274 ENDIF 275 276 277 !----------------------------------------------------------------------- 278 ! 5.calcul des traansmissions depuis le sol, cas diffus: 279 ! ------------------------------------------------------ 280 281 DO l=1,nlevel 282 DO ig=1,ncount 283 ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66) 284 ENDDO 285 ENDDO 286 287 IF (lwrite) THEN 288 PRINT* 289 PRINT*,'Diagnostique des taux de chauffage solaires' 290 PRINT*,' 3 transmission avec les sol' 291 PRINT*,'niveau transmission' 292 DO l=1,nlevel 293 PRINT*,l,ztrref(igout,l) 294 ENDDO 295 ENDIF 296 297 !----------------------------------------------------------------------- 298 ! 6.ajout a l'echauffement de la contribution du ray. sol. reflechit: 299 ! ------------------------------------------------------------------- 300 301 DO l=1,nlayer 302 DO ig=1,ncount 303 zdtsw(ig,l)=zdtsw(ig,l)+ & 304 g*zfsrfref(ig)*(ztrref(ig,l+1)-ztrref(ig,l))/ & 305 (cpp*(zplev(ig,l+1)-zplev(ig,l))) 306 ENDDO 307 ENDDO 308 309 !----------------------------------------------------------------------- 310 ! 10. sorites eventuelles: 102 311 ! ------------------------ 103 312 104 alph=phaut/prad 105 DO j=1,npts 106 ! !!!!!! 107 pmu(j)=sqrt(1224.*pmu(j)*pmu(j)+1.)/35. 108 ! $ (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j)) 109 END DO 110 111 END SUBROUTINE mucorr 112 113 END MODULE radiative 313 IF (lwrite) THEN 314 PRINT* 315 PRINT*,'Diagnostique des taux de chauffage solaires:' 316 PRINT*,' 3 taux de chauffage total' 317 DO l=1,nlayer 318 PRINT*,zdtsw(igout,l) 319 ENDDO 320 ENDIF 321 322 IF (ldiurn) THEN 323 fsrfvis(:)=0. 324 CALL monscatter(ncount,fsrfvis,index,zflux) 325 dtsw(:,:)=0. 326 DO l=1,nlayer 327 CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l)) 328 ENDDO 329 ELSE 330 print*,'NOT DIURNE' 331 fsrfvis(:)=zflux(:) 332 dtsw(:,:)=zdtsw(:,:) 333 ENDIF 334 ! call dump2d(iim,jjm-1,zflux(2),'ZFLUX ') 335 ! call dump2d(iim,jjm-1,fsrfvis(2),'FSRVIS ') 336 ! call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') 337 ! call dump2d(iim,jjm-1,pmu(2),'pmu ') 338 ! call dump2d(iim,jjm-1,pfract(2),'pfract ') 339 ! call dump2d(iim,jjm-1,albedo(2),'albedo ') 340 ! call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') 341 342 343 END SUBROUTINE sw 344 345 END MODULE radiative
Note: See TracChangeset
for help on using the changeset viewer.
