Changeset 4201 for dynamico_lmdz/simple_physics/phyparam
- Timestamp:
- Dec 20, 2019, 2:25:23 PM (5 years ago)
- Location:
- dynamico_lmdz/simple_physics/phyparam
- Files:
-
- 1 deleted
- 2 edited
- 1 moved
-
param/phyparam.F (modified) (2 diffs)
-
param/zenang.F (deleted)
-
physics/radiative_sw.F90 (moved) (moved from dynamico_lmdz/simple_physics/phyparam/physics/radiative.F90) (3 diffs)
-
physics/solar.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
dynamico_lmdz/simple_physics/phyparam/param/phyparam.F
r4200 r4201 14 14 USE astronomy 15 15 USE vdif_mod, ONLY : vdif 16 USE solar, ONLY : solang 17 USE radiative , ONLY : mucorr,sw16 USE solar, ONLY : solang, zenang, mucorr 17 USE radiative_sw, ONLY : sw 18 18 USE radiative_lw, ONLY : lw 19 19 … … 373 373 else 374 374 zdtime=ptimestep*float(iradia) 375 CALL zenang( zls,gmtime,zdtime,lati,long,mu0,fract)375 CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract) 376 376 print*,'ZENANG ' 377 377 endif -
dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90
r4199 r4201 1 MODULE radiative 1 MODULE radiative_sw 2 2 IMPLICIT NONE 3 3 SAVE … … 5 5 PRIVATE 6 6 7 PUBLIC :: mucorr, sw, lwtr7 PUBLIC :: sw 8 8 9 9 CONTAINS 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 whithout16 ! diurnal cycle.17 !18 ! parmeters :19 ! -----------20 !21 ! Input :22 ! -------23 ! npts number of points24 ! pdeclin solar declinaison25 ! plat(npts) latitude26 ! phaut hauteur typique de l'atmosphere27 ! prad rayon planetaire28 !29 ! Output :30 ! --------31 ! pmu(npts) equivalent cosinus of the solar angle32 ! pfract(npts) fractionnal day33 !34 !=======================================================================35 36 !-----------------------------------------------------------------------37 !38 ! 0. Declarations :39 ! -----------------40 41 ! Arguments :42 ! -----------43 INTEGER, INTENT(IN) :: npts44 REAL, INTENT(IN) :: phaut, prad, pdeclin, plat(npts)45 REAL, INTENT(OUT) :: pmu(npts), pfract(npts)46 !47 ! Local variables :48 ! -----------------49 INTEGER j50 REAL z,cz,sz,tz,phi,cphi,sphi,tphi51 REAL ap,a,t,b52 REAL alph53 54 REAL, PARAMETER :: pi=2.*asin(1.)55 56 !-----------------------------------------------------------------------57 58 ! print*,'npts,pdeclin',npts,pdeclin59 z = pdeclin60 cz = cos (z)61 sz = sin (z)62 ! print*,'cz,sz',cz,sz63 64 DO j = 1, npts65 66 phi = plat(j)67 cphi = cos(phi)68 if (cphi.le.1.e-9) cphi=1.e-969 sphi = sin(phi)70 tphi = sphi / cphi71 b = cphi * cz72 t = -tphi * sz / cz73 a = 1.0 - t*t74 ap = a75 76 IF(t.eq.0.) then77 t=0.5*pi78 ELSE79 IF (a.lt.0.) a = 0.80 t = sqrt(a) / t81 IF (t.lt.0.) then82 t = -atan (-t) + pi83 ELSE84 t = atan(t)85 ENDIF86 ENDIF87 88 pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi89 pfract(j) = t / pi90 IF (ap .lt.0.) then91 pmu(j) = sphi * sz92 pfract(j) = 1.093 ENDIF94 IF (pmu(j).le.0.0) pmu(j) = 0.095 pmu(j) = pmu(j) / pfract(j)96 IF (pmu(j).eq.0.) pfract(j) = 0.97 98 END DO99 100 !-----------------------------------------------------------------------101 ! correction de rotondite:102 ! ------------------------103 104 alph=phaut/prad105 DO j=1,npts106 ! !!!!!!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 DO110 111 END SUBROUTINE mucorr112 10 113 11 PURE SUBROUTINE monGATHER(n,a,b,index) … … 128 26 lwrite) 129 27 USE phys_const 130 !======================================================================= 131 ! 132 ! Rayonnement solaire en atmosphere non diffusante avec un 133 ! coefficient d'absoprption gris. 134 ! 135 !======================================================================= 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 187 ENDIF 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: 311 ! ------------------------ 312 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 PURE SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm) 346 INTEGER, INTENT(IN) :: ngrid 347 REAL, INTENT(IN) :: coef 348 LOGICAL, INTENT(IN) :: lstrong 349 REAL, INTENT(IN) :: dup(ngrid) 350 REAL, INTENT(OUT) :: transm(ngrid) 351 INTEGER ig 352 IF(lstrong) THEN 353 DO ig=1,ngrid 354 transm(ig)=exp(-coef*sqrt(dup(ig))) 355 ENDDO 356 ELSE 357 DO ig=1,ngrid 358 transm(ig)=exp(-coef*dup(ig)) 359 ENDDO 360 ENDIF 361 362 END SUBROUTINE lwtr 363 364 END MODULE radiative 28 !======================================================================= 29 ! 30 ! Rayonnement solaire en atmosphere non diffusante avec un 31 ! coefficient d'absoprption gris. 32 ! 33 !======================================================================= 34 ! 35 ! declarations: 36 ! ------------- 37 ! 38 ! 39 ! arguments: 40 ! ---------- 41 ! 42 INTEGER ngrid,nlayer 43 REAL albedo(ngrid),coefvis 44 REAL pmu(ngrid),pfract(ngrid) 45 REAL plevel(ngrid,nlayer+1),ps_rad 46 REAL psolarf0 47 REAL fsrfvis(ngrid),dtsw(ngrid,nlayer) 48 LOGICAL lwrite,ldiurn 49 ! 50 ! variables locales: 51 ! ------------------ 52 ! 53 54 REAL zalb(ngrid),zmu(ngrid),zfract(ngrid) 55 REAL zplev(ngrid,nlayer+1) 56 REAL zflux(ngrid),zdtsw(ngrid,nlayer) 57 58 INTEGER ig,l,nlevel,index(ngrid),ncount,igout 59 REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1) 60 REAL zfsrfref(ngrid) 61 REAL z1(ngrid) 62 REAL zu(ngrid,nlayer+1) 63 REAL tau0 64 65 !----------------------------------------------------------------------- 66 ! 1. initialisations: 67 ! ------------------- 68 69 70 nlevel=nlayer+1 71 72 !----------------------------------------------------------------------- 73 ! Definitions des tableaux locaux pour les points ensoleilles: 74 ! ------------------------------------------------------------ 75 76 IF (ldiurn) THEN 77 ncount=0 78 DO ig=1,ngrid 79 index(ig)=0 80 ENDDO 81 DO ig=1,ngrid 82 IF(pfract(ig).GT.1.e-6) THEN 83 ncount=ncount+1 84 index(ncount)=ig 85 ENDIF 86 ENDDO 87 CALL monGATHER(ncount,zfract,pfract,index) 88 CALL monGATHER(ncount,zmu,pmu,index) 89 CALL monGATHER(ncount,zalb,albedo,index) 90 DO l=1,nlevel 91 CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index) 92 ENDDO 93 ELSE 94 ncount=ngrid 95 zfract(:)=pfract(:) 96 zmu(:)=pmu(:) 97 zalb(:)=albedo(:) 98 zplev(:,:)=plevel(:,:) 99 ENDIF 100 101 !----------------------------------------------------------------------- 102 ! calcul des profondeurs optiques integres depuis p=0: 103 ! ---------------------------------------------------- 104 105 tau0=-.5*log(coefvis) 106 107 ! calcul de la partie homogene de l'opacite 108 tau0=tau0/ps_rad 109 DO l=1,nlayer+1 110 DO ig=1,ncount 111 zu(ig,l)=tau0*zplev(ig,l) 112 ENDDO 113 ENDDO 114 115 !----------------------------------------------------------------------- 116 ! 2. calcul de la transmission depuis le sommet de l'atmosphere: 117 ! ----------------------------------------------------------- 118 119 DO l=1,nlevel 120 DO ig=1,ncount 121 ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig)) 122 ENDDO 123 ENDDO 124 125 IF (lwrite) THEN 126 igout=ncount/2+1 127 PRINT* 128 PRINT*,'Diagnostique des transmission dans le spectre solaire' 129 PRINT*,'zfract, zmu, zalb' 130 PRINT*,zfract(igout), zmu(igout), zalb(igout) 131 PRINT*,'Pression, quantite d abs, transmission' 132 DO l=1,nlayer+1 133 PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l) 134 ENDDO 135 ENDIF 136 137 !----------------------------------------------------------------------- 138 ! 3. taux de chauffage, ray. solaire direct: 139 ! ------------------------------------------ 140 141 DO l=1,nlayer 142 DO ig=1,ncount 143 zdtsw(ig,l)=g*psolarf0*zfract(ig)*zmu(ig)* & 144 (ztrdir(ig,l+1)-ztrdir(ig,l))/ & 145 (cpp*(zplev(ig,l)-zplev(ig,l+1))) 146 ENDDO 147 ENDDO 148 IF (lwrite) THEN 149 PRINT* 150 PRINT*,'Diagnostique des taux de chauffage solaires:' 151 PRINT*,' 1 taux de chauffage lie au ray. solaire direct' 152 DO l=1,nlayer 153 PRINT*,zdtsw(igout,l) 154 ENDDO 155 ENDIF 156 157 158 !----------------------------------------------------------------------- 159 ! 4. calcul du flux solaire arrivant sur le sol: 160 ! ---------------------------------------------- 161 162 DO ig=1,ncount 163 z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1) 164 zflux(ig)=(1.-zalb(ig))*z1(ig) 165 zfsrfref(ig)= zalb(ig)*z1(ig) 166 ENDDO 167 IF (lwrite) THEN 168 PRINT* 169 PRINT*,'Diagnostique des taux de chauffage solaires:' 170 PRINT*,' 2 flux solaire net incident sur le sol' 171 PRINT*,zflux(igout) 172 ENDIF 173 174 175 !----------------------------------------------------------------------- 176 ! 5.calcul des traansmissions depuis le sol, cas diffus: 177 ! ------------------------------------------------------ 178 179 DO l=1,nlevel 180 DO ig=1,ncount 181 ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66) 182 ENDDO 183 ENDDO 184 185 IF (lwrite) THEN 186 PRINT* 187 PRINT*,'Diagnostique des taux de chauffage solaires' 188 PRINT*,' 3 transmission avec les sol' 189 PRINT*,'niveau transmission' 190 DO l=1,nlevel 191 PRINT*,l,ztrref(igout,l) 192 ENDDO 193 ENDIF 194 195 !----------------------------------------------------------------------- 196 ! 6.ajout a l'echauffement de la contribution du ray. sol. reflechit: 197 ! ------------------------------------------------------------------- 198 199 DO l=1,nlayer 200 DO ig=1,ncount 201 zdtsw(ig,l)=zdtsw(ig,l)+ & 202 g*zfsrfref(ig)*(ztrref(ig,l+1)-ztrref(ig,l))/ & 203 (cpp*(zplev(ig,l+1)-zplev(ig,l))) 204 ENDDO 205 ENDDO 206 207 !----------------------------------------------------------------------- 208 ! 10. sorites eventuelles: 209 ! ------------------------ 210 211 IF (lwrite) THEN 212 PRINT* 213 PRINT*,'Diagnostique des taux de chauffage solaires:' 214 PRINT*,' 3 taux de chauffage total' 215 DO l=1,nlayer 216 PRINT*,zdtsw(igout,l) 217 ENDDO 218 ENDIF 219 220 IF (ldiurn) THEN 221 fsrfvis(:)=0. 222 CALL monscatter(ncount,fsrfvis,index,zflux) 223 dtsw(:,:)=0. 224 DO l=1,nlayer 225 CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l)) 226 ENDDO 227 ELSE 228 print*,'NOT DIURNE' 229 fsrfvis(:)=zflux(:) 230 dtsw(:,:)=zdtsw(:,:) 231 ENDIF 232 ! call dump2d(iim,jjm-1,zflux(2),'ZFLUX ') 233 ! call dump2d(iim,jjm-1,fsrfvis(2),'FSRVIS ') 234 ! call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') 235 ! call dump2d(iim,jjm-1,pmu(2),'pmu ') 236 ! call dump2d(iim,jjm-1,pfract(2),'pfract ') 237 ! call dump2d(iim,jjm-1,albedo(2),'albedo ') 238 ! call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') 239 240 241 END SUBROUTINE sw 242 243 END MODULE radiative_sw -
dynamico_lmdz/simple_physics/phyparam/physics/solar.F90
r4200 r4201 2 2 IMPLICIT NONE 3 3 4 PRIVATE 5 6 REAL, PARAMETER :: pi_local = 4.0 * ATAN(1.0), deux_pi_local = 2.0 * pi_local 7 8 PUBLIC :: solang, zenang, mucorr 9 4 10 CONTAINS 5 11 6 subroutinesolang ( kgrid,psilon,pcolon,psilat,pcolat, &12 PURE SUBROUTINE solang ( kgrid,psilon,pcolon,psilat,pcolat, & 7 13 ptim1,ptim2,ptim3,pmu0,pfract ) 8 9 ! 10 !**** *LW* - ORGANIZES THE LONGWAVE CALCULATIONS 11 ! 12 ! PURPOSE. 13 ! -------- 14 ! CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID 15 ! 16 !** INTERFACE. 17 ! ---------- 18 ! SUBROUTINE SOLANG ( KGRID ) 19 ! 20 ! EXPLICIT ARGUMENTS : 21 ! -------------------- 22 ! ==== INPUTS === 23 ! 24 ! PSILON(KGRID) : SINUS OF THE LONGITUDE 25 ! PCOLON(KGRID) : COSINUS OF THE LONGITUDE 26 ! PSILAT(KGRID) : SINUS OF THE LATITUDE 27 ! PCOLAT(KGRID) : COSINUS OF THE LATITUDE 28 ! PTIM1 : SIN(DECLI) 29 ! PTIM2 : COS(DECLI)*COS(TIME) 30 ! PTIM3 : SIN(DECLI)*SIN(TIME) 31 ! 32 ! ==== OUTPUTS === 33 ! 34 ! PMU0 (KGRID) : SOLAR ANGLE 35 ! PFRACT(KGRID) : DAY FRACTION OF THE TIME INTERVAL 36 ! 37 ! IMPLICIT ARGUMENTS : NONE 38 ! -------------------- 39 ! 40 ! METHOD. 41 ! ------- 42 ! 43 ! EXTERNALS. 44 ! ---------- 45 ! 46 ! NONE 47 ! 48 ! REFERENCE. 49 ! ---------- 50 ! 51 ! RADIATIVE PROCESSES IN METEOROLOGIE AND CLIMATOLOGIE 52 ! PALTRIDGE AND PLATT 53 ! 54 ! AUTHOR. 55 ! ------- 56 ! FREDERIC HOURDIN 57 ! 58 ! MODIFICATIONS. 59 ! -------------- 60 ! ORIGINAL :90-01-14 61 ! 92-02-14 CALCULATIONS DONE THE ENTIER GRID (J.Polcher) 62 !----------------------------------------------------------------------- 63 ! 64 ! ------------------------------------------------------------------ 65 !----------------------------------------------------------------------- 66 ! 67 !* 0.1 ARGUMENTS 68 ! --------- 69 ! 70 INTEGER kgrid 71 REAL ptim1,ptim2,ptim3 72 REAL psilon(kgrid),pcolon(kgrid),pmu0(kgrid),pfract(kgrid) 73 REAL psilat(kgrid), pcolat(kgrid) 74 75 INTEGER jl 76 REAL ztim1,ztim2,ztim3 77 !------------------------------------------------------------------------ 78 !------------------------------------------------------------------------ 79 !------------------------------------------------------------------------ 80 ! 81 !------------------------------------------------------------------------ 82 ! 83 !* 1. INITIALISATION 84 ! -------------- 85 ! 86 100 CONTINUE 87 ! 88 DO jl=1,kgrid 89 pmu0(jl)=0. 90 pfract(jl)=0. 91 ENDDO 92 ! 93 !* 1.1 COMPUTATION OF THE SOLAR ANGLE 94 ! ------------------------------ 95 ! 96 DO jl=1,kgrid 97 ztim1=psilat(jl)*ptim1 98 ztim2=pcolat(jl)*ptim2 99 ztim3=pcolat(jl)*ptim3 100 pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl) 101 ENDDO 102 ! 103 !* 1.2 DISTINCTION BETWEEN DAY AND NIGHT 104 ! --------------------------------- 105 ! 106 DO jl=1,kgrid 107 IF (pmu0(jl).gt.0.) THEN 14 15 !----------------------------------------------------------------------- 16 ! CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID 17 ! 18 ! ==== INPUTS === 19 ! 20 ! PSILON(KGRID) : SINUS OF THE LONGITUDE 21 ! PCOLON(KGRID) : COSINUS OF THE LONGITUDE 22 ! PSILAT(KGRID) : SINUS OF THE LATITUDE 23 ! PCOLAT(KGRID) : COSINUS OF THE LATITUDE 24 ! PTIM1 : SIN(DECLI) 25 ! PTIM2 : COS(DECLI)*COS(TIME) 26 ! PTIM3 : SIN(DECLI)*SIN(TIME) 27 ! 28 ! ==== OUTPUTS === 29 ! 30 ! PMU0 (KGRID) : SOLAR ANGLE 31 ! PFRACT(KGRID) : DAY FRACTION OF THE TIME INTERVAL 32 ! 33 ! REFERENCE. 34 ! ---------- 35 ! 36 ! RADIATIVE PROCESSES IN METEOROLOGIE AND CLIMATOLOGIE 37 ! PALTRIDGE AND PLATT 38 ! 39 ! AUTHOR. 40 ! ------- 41 ! FREDERIC HOURDIN 42 ! 43 ! MODIFICATIONS. 44 ! -------------- 45 ! ORIGINAL :90-01-14 46 ! 92-02-14 CALCULATIONS DONE THE ENTIER GRID (J.Polcher) 47 !----------------------------------------------------------------------- 48 49 INTEGER, INTENT(IN) :: kgrid 50 REAL, INTENT(IN) :: ptim1,ptim2,ptim3 51 REAL, INTENT(IN) :: psilon(kgrid), pcolon(kgrid) 52 REAL, INTENT(IN) :: psilat(kgrid), pcolat(kgrid) 53 REAL, INTENT(OUT) :: pmu0(kgrid), pfract(kgrid) 54 55 INTEGER jl 56 REAL ztim1,ztim2,ztim3 57 !------------------------------------------------------------------------ 58 !------------------------------------------------------------------------ 59 !------------------------------------------------------------------------ 60 ! 61 !------------------------------------------------------------------------ 62 ! 63 !* 1. INITIALISATION 64 ! -------------- 65 ! 66 ! 67 DO jl=1,kgrid 68 pmu0(jl)=0. 69 pfract(jl)=0. 70 ENDDO 71 ! 72 !* 1.1 COMPUTATION OF THE SOLAR ANGLE 73 ! ------------------------------ 74 ! 75 DO jl=1,kgrid 76 ztim1=psilat(jl)*ptim1 77 ztim2=pcolat(jl)*ptim2 78 ztim3=pcolat(jl)*ptim3 79 pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl) 80 ENDDO 81 ! 82 !* 1.2 DISTINCTION BETWEEN DAY AND NIGHT 83 ! --------------------------------- 84 ! 85 DO jl=1,kgrid 86 IF (pmu0(jl).gt.0.) THEN 108 87 pfract(jl)=1. 109 88 ELSE 110 89 pmu0(jl)=0. 111 pfract(jl)=0. 112 ENDIF 113 ENDDO 114 ! 115 116 END subroutine solang 117 END MODULE solar 90 pfract(jl)=0. 91 ENDIF 92 ENDDO 93 94 END SUBROUTINE solang 95 96 SUBROUTINE zenang(klon,longi,gmtime,pdtrad,lat,long, & 97 & pmu0,frac) 98 USE astronomy 99 100 !============================================================= 101 ! Auteur : O. Boucher (LMD/CNRS) 102 ! d'apres les routines zenith et angle de Z.X. Li 103 ! Objet : calculer les valeurs moyennes du cos de l'angle zenithal 104 ! et l'ensoleillement moyen entre gmtime1 et gmtime2 105 ! connaissant la declinaison, la latitude et la longitude. 106 ! Rque : Different de la routine angle en ce sens que zenang 107 ! fournit des moyennes de pmu0 et non des valeurs 108 ! instantanees, du coup frac prend toutes les valeurs 109 ! entre 0 et 1. 110 ! Date : premiere version le 13 decembre 1994 111 ! revu pour GCM le 30 septembre 1996 112 !=============================================================== 113 ! longi----INPUT : la longitude vraie de la terre dans son plan 114 ! solaire a partir de l'equinoxe de printemps (degre) 115 ! gmtime---INPUT : temps universel en fraction de jour 116 ! pdtrad---INPUT : pas de temps du rayonnement (secondes) 117 ! lat------INPUT : latitude en degres 118 ! long-----INPUT : longitude en degres 119 ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad 120 ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad 121 !================================================================ 122 integer klon 123 !================================================================ 124 real longi, gmtime, pdtrad 125 real lat(klon), long(klon), pmu0(klon), frac(klon) 126 !================================================================ 127 integer i 128 real gmtime1, gmtime2 129 real incl 130 real omega1, omega2, omega 131 ! omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi. 132 ! omega : heure en radian du coucher de soleil 133 ! -omega est donc l'heure en radian de lever du soleil 134 real omegadeb, omegafin 135 real zfrac1, zfrac2, z1_mu, z2_mu 136 ! declinaison en radian 137 real lat_sun 138 ! longitude solaire en radian 139 real lon_sun 140 ! latitude du pt de grille en radian 141 real latr 142 !================================================================ 143 ! 144 ! incl=R_incl * pi_local / 180. 145 print*,'Obliquite =' ,obliquit 146 incl=obliquit * pi_local / 180. 147 ! 148 ! lon_sun = longi * pi_local / 180.0 149 lon_sun = longi 150 lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) 151 ! 152 gmtime1=gmtime*86400. 153 gmtime2=gmtime*86400.+pdtrad 154 ! 155 DO i = 1, klon 156 ! 157 ! latr = lat(i) * pi_local / 180. 158 latr = lat(i) 159 ! 160 !--pose probleme quand lat=+/-90 degres 161 ! 162 ! omega = -TAN(latr)*TAN(lat_sun) 163 ! omega = ACOS(omega) 164 ! IF (latr.GE.(pi_local/2.+lat_sun) 165 ! . .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN 166 ! omega = 0.0 ! nuit polaire 167 ! ENDIF 168 ! IF (latr.GE.(pi_local/2.-lat_sun) 169 ! . .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN 170 ! omega = pi_local ! journee polaire 171 ! ENDIF 172 ! 173 !--remplace par cela (le cas par defaut est different) 174 ! 175 !--nuit polaire 176 omega=0.0 177 IF (latr.GE.(pi_local/2.-lat_sun) & 178 & .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN 179 ! journee polaire 180 omega = pi_local 181 ENDIF 182 IF (latr.LT.(pi_local/2.+lat_sun).AND. & 183 & latr.GT.(-pi_local/2.+lat_sun).AND. & 184 & latr.LT.(pi_local/2.-lat_sun).AND. & 185 & latr.GT.(-pi_local/2.-lat_sun)) THEN 186 omega = -TAN(latr)*TAN(lat_sun) 187 omega = ACOS(omega) 188 ENDIF 189 ! 190 omega1 = gmtime1 + long(i)*86400.0/360.0 191 omega1 = omega1 / 86400.0*deux_pi_local 192 omega1 = MOD (omega1+deux_pi_local, deux_pi_local) 193 omega1 = omega1 - pi_local 194 ! 195 omega2 = gmtime2 + long(i)*86400.0/360.0 196 omega2 = omega2 / 86400.0*deux_pi_local 197 omega2 = MOD (omega2+deux_pi_local, deux_pi_local) 198 omega2 = omega2 - pi_local 199 ! 200 !--on est dans la meme journee locale 201 IF (omega1.LE.omega2) THEN 202 ! 203 IF (omega2.LE.-omega .OR. omega1.GE.omega .OR. omega.LT.1e-5) & 204 THEN 205 !--nuit 206 frac(i)=0.0 207 pmu0(i)=0.0 208 !--jour+nuit/jou 209 ELSE 210 omegadeb=MAX(-omega,omega1) 211 omegafin=MIN(omega,omega2) 212 frac(i)=(omegafin-omegadeb)/(omega2-omega1) 213 pmu0(i)=SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)* & 214 (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb) 215 ENDIF 216 ! 217 !---omega1 GT omega2 -- a cheval sur deux journees 218 ELSE 219 ! 220 !-------------------entre omega1 et pi 221 !--nuit 222 IF (omega1.GE.omega) THEN 223 zfrac1=0.0 224 z1_mu =0.0 225 !--jour+nuit 226 ELSE 227 omegadeb=MAX(-omega,omega1) 228 omegafin=omega 229 zfrac1=omegafin-omegadeb 230 z1_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)* & 231 (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb) 232 ENDIF 233 !---------------------entre -pi et omega2 234 !--nuit 235 IF (omega2.LE.-omega) THEN 236 zfrac2=0.0 237 z2_mu =0.0 238 !--jour+nuit 239 ELSE 240 omegadeb=-omega 241 omegafin=MIN(omega,omega2) 242 zfrac2=omegafin-omegadeb 243 z2_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)* & 244 (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb) 245 ! 246 ENDIF 247 !-----------------------moyenne 248 frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1) 249 pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10) 250 ! 251 !---comparaison omega1 et omega2 252 ENDIF 253 ! 254 ENDDO 255 ! 256 END SUBROUTINE zenang 257 258 PURE SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad) 259 260 !======================================================================= 261 ! 262 ! Calcul of equivalent solar angle and and fraction of day whithout 263 ! diurnal cycle. 264 ! 265 ! parmeters : 266 ! ----------- 267 ! 268 ! Input : 269 ! ------- 270 ! npts number of points 271 ! pdeclin solar declinaison 272 ! plat(npts) latitude 273 ! phaut hauteur typique de l'atmosphere 274 ! prad rayon planetaire 275 ! 276 ! Output : 277 ! -------- 278 ! pmu(npts) equivalent cosinus of the solar angle 279 ! pfract(npts) fractionnal day 280 ! 281 !======================================================================= 282 283 !----------------------------------------------------------------------- 284 ! 285 ! 0. Declarations : 286 ! ----------------- 287 288 ! Arguments : 289 ! ----------- 290 INTEGER, INTENT(IN) :: npts 291 REAL, INTENT(IN) :: phaut, prad, pdeclin, plat(npts) 292 REAL, INTENT(OUT) :: pmu(npts), pfract(npts) 293 ! 294 ! Local variables : 295 ! ----------------- 296 INTEGER j 297 REAL z,cz,sz,tz,phi,cphi,sphi,tphi 298 REAL ap,a,t,b 299 REAL alph 300 301 REAL, PARAMETER :: pi=2.*asin(1.) 302 303 !----------------------------------------------------------------------- 304 305 ! print*,'npts,pdeclin',npts,pdeclin 306 z = pdeclin 307 cz = cos (z) 308 sz = sin (z) 309 ! print*,'cz,sz',cz,sz 310 311 DO j = 1, npts 312 313 phi = plat(j) 314 cphi = cos(phi) 315 if (cphi.le.1.e-9) cphi=1.e-9 316 sphi = sin(phi) 317 tphi = sphi / cphi 318 b = cphi * cz 319 t = -tphi * sz / cz 320 a = 1.0 - t*t 321 ap = a 322 323 IF(t.eq.0.) then 324 t=0.5*pi 325 ELSE 326 IF (a.lt.0.) a = 0. 327 t = sqrt(a) / t 328 IF (t.lt.0.) then 329 t = -atan (-t) + pi 330 ELSE 331 t = atan(t) 332 ENDIF 333 ENDIF 334 335 pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi 336 pfract(j) = t / pi 337 IF (ap .lt.0.) then 338 pmu(j) = sphi * sz 339 pfract(j) = 1.0 340 ENDIF 341 IF (pmu(j).le.0.0) pmu(j) = 0.0 342 pmu(j) = pmu(j) / pfract(j) 343 IF (pmu(j).eq.0.) pfract(j) = 0. 344 345 END DO 346 347 !----------------------------------------------------------------------- 348 ! correction de rotondite: 349 ! ------------------------ 350 351 alph=phaut/prad 352 DO j=1,npts 353 ! !!!!!! 354 pmu(j)=sqrt(1224.*pmu(j)*pmu(j)+1.)/35. 355 ! $ (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j)) 356 END DO 357 358 END SUBROUTINE mucorr 359 360 END MODULE solar
Note: See TracChangeset
for help on using the changeset viewer.
