Index: namico_lmdz/simple_physics/phyparam/param/gather.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/gather.F (revision 4197) +++ (revision ) @@ -1,16 +1,0 @@ - SUBROUTINE monGATHER(n,a,b,index) -c - IMPLICIT NONE -C - INTEGER n,index(n),i - REAL a(n),b(n) -c - DO 100 i=1,n - a(i)=b(index(i)) -100 CONTINUE -c - RETURN - END -c -c voir aussi vec_gather(v,vindices,count,r)...p.11-14 -c Index: /dynamico_lmdz/simple_physics/phyparam/param/phyparam.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/phyparam.F (revision 4197) +++ /dynamico_lmdz/simple_physics/phyparam/param/phyparam.F (revision 4198) @@ -14,5 +14,5 @@ USE astronomy USE vdif_mod, ONLY : vdif - USE radiative, ONLY : mucorr + USE radiative, ONLY : mucorr, sw c IMPLICIT NONE Index: namico_lmdz/simple_physics/phyparam/param/sw.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/sw.F (revision 4197) +++ (revision ) @@ -1,227 +1,0 @@ - SUBROUTINE sw(ngrid,nlayer,ldiurn, - $ coefvis,albedo, - $ plevel,ps_rad,pmu,pfract,psolarf0, - $ fsrfvis,dtsw, - $ lwrite) - USE phys_const - IMPLICIT NONE -c======================================================================= -c -c Rayonnement solaire en atmosphere non diffusante avec un -c coefficient d'absoprption gris. -c -c======================================================================= -c -c declarations: -c ------------- -c -c -c arguments: -c ---------- -c - INTEGER ngrid,nlayer - REAL albedo(ngrid),coefvis - REAL pmu(ngrid),pfract(ngrid) - REAL plevel(ngrid,nlayer+1),ps_rad - REAL psolarf0 - REAL fsrfvis(ngrid),dtsw(ngrid,nlayer) - LOGICAL lwrite,ldiurn -c -c variables locales: -c ------------------ -c - - REAL zalb(ngrid),zmu(ngrid),zfract(ngrid) - REAL zplev(ngrid,nlayer+1) - REAL zflux(ngrid),zdtsw(ngrid,nlayer) - - INTEGER ig,l,nlevel,index(ngrid),ncount,igout - REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1) - REAL zfsrfref(ngrid) - REAL z1(ngrid) - REAL zu(ngrid,nlayer+1) - REAL tau0 - - LOGICAL firstcall - SAVE firstcall - DATA firstcall/.true./ -!$OMP THREADPRIVATE(firstcall) - -c----------------------------------------------------------------------- -c 1. initialisations: -c ------------------- - - - nlevel=nlayer+1 - -c----------------------------------------------------------------------- -c Definitions des tableaux locaux pour les points ensoleilles: -c ------------------------------------------------------------ - - IF (ldiurn) THEN - ncount=0 - DO ig=1,ngrid - index(ig)=0 - ENDDO - DO ig=1,ngrid - IF(pfract(ig).GT.1.e-6) THEN - ncount=ncount+1 - index(ncount)=ig - ENDIF - ENDDO - CALL monGATHER(ncount,zfract,pfract,index) - CALL monGATHER(ncount,zmu,pmu,index) - CALL monGATHER(ncount,zalb,albedo,index) - DO l=1,nlevel - CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index) - ENDDO - ELSE - ncount=ngrid - zfract(:)=pfract(:) - zmu(:)=pmu(:) - zalb(:)=albedo(:) - zplev(:,:)=plevel(:,:) - ENDIF - -c----------------------------------------------------------------------- -c calcul des profondeurs optiques integres depuis p=0: -c ---------------------------------------------------- - - tau0=-.5*log(coefvis) - -c calcul de la partie homogene de l'opacite - tau0=tau0/ps_rad - DO l=1,nlayer+1 - DO ig=1,ncount - zu(ig,l)=tau0*zplev(ig,l) - ENDDO - ENDDO - -c----------------------------------------------------------------------- -c 2. calcul de la transmission depuis le sommet de l'atmosphere: -c ----------------------------------------------------------- - - DO l=1,nlevel - DO ig=1,ncount - ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig)) - ENDDO - ENDDO - - IF (lwrite) THEN - igout=ncount/2+1 - PRINT* - PRINT*,'Diagnostique des transmission dans le spectre solaire' - PRINT*,'zfract, zmu, zalb' - PRINT*,zfract(igout), zmu(igout), zalb(igout) - PRINT*,'Pression, quantite d abs, transmission' - DO l=1,nlayer+1 - PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l) - ENDDO - ENDIF - -c----------------------------------------------------------------------- -c 3. taux de chauffage, ray. solaire direct: -c ------------------------------------------ - - DO l=1,nlayer - DO ig=1,ncount - zdtsw(ig,l)=g*psolarf0*zfract(ig)*zmu(ig)* - $ (ztrdir(ig,l+1)-ztrdir(ig,l))/ - $ (cpp*(zplev(ig,l)-zplev(ig,l+1))) - ENDDO - ENDDO - IF (lwrite) THEN - PRINT* - PRINT*,'Diagnostique des taux de chauffage solaires:' - PRINT*,' 1 taux de chauffage lie au ray. solaire direct' - DO l=1,nlayer - PRINT*,zdtsw(igout,l) - ENDDO - ENDIF - - -c----------------------------------------------------------------------- -c 4. calcul du flux solaire arrivant sur le sol: -c ---------------------------------------------- - - DO ig=1,ncount - z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1) - zflux(ig)=(1.-zalb(ig))*z1(ig) - zfsrfref(ig)= zalb(ig)*z1(ig) - ENDDO - IF (lwrite) THEN - PRINT* - PRINT*,'Diagnostique des taux de chauffage solaires:' - PRINT*,' 2 flux solaire net incident sur le sol' - PRINT*,zflux(igout) - ENDIF - - -c----------------------------------------------------------------------- -c 5.calcul des traansmissions depuis le sol, cas diffus: -c ------------------------------------------------------ - - DO l=1,nlevel - DO ig=1,ncount - ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66) - ENDDO - ENDDO - - IF (lwrite) THEN - PRINT* - PRINT*,'Diagnostique des taux de chauffage solaires' - PRINT*,' 3 transmission avec les sol' - PRINT*,'niveau transmission' - DO l=1,nlevel - PRINT*,l,ztrref(igout,l) - ENDDO - ENDIF - -c----------------------------------------------------------------------- -c 6.ajout a l'echauffement de la contribution du ray. sol. reflechit: -c ------------------------------------------------------------------- - - DO l=1,nlayer - DO ig=1,ncount - zdtsw(ig,l)=zdtsw(ig,l)+ - $ g*zfsrfref(ig)*(ztrref(ig,l+1)-ztrref(ig,l))/ - $ (cpp*(zplev(ig,l+1)-zplev(ig,l))) - ENDDO - ENDDO - -c----------------------------------------------------------------------- -c 10. sorites eventuelles: -c ------------------------ - - IF (lwrite) THEN - PRINT* - PRINT*,'Diagnostique des taux de chauffage solaires:' - PRINT*,' 3 taux de chauffage total' - DO l=1,nlayer - PRINT*,zdtsw(igout,l) - ENDDO - ENDIF - - IF (ldiurn) THEN - fsrfvis(:)=0. - CALL monscatter(ncount,fsrfvis,index,zflux) - dtsw(:,:)=0. - DO l=1,nlayer - CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l)) - ENDDO - ELSE - print*,'NOT DIURNE' - fsrfvis(:)=zflux(:) - dtsw(:,:)=zdtsw(:,:) - ENDIF -c call dump2d(iim,jjm-1,zflux(2),'ZFLUX ') -c call dump2d(iim,jjm-1,fsrfvis(2),'FSRVIS ') -c call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') -c call dump2d(iim,jjm-1,pmu(2),'pmu ') -c call dump2d(iim,jjm-1,pfract(2),'pfract ') -c call dump2d(iim,jjm-1,albedo(2),'albedo ') -c call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') - - - RETURN - END Index: /dynamico_lmdz/simple_physics/phyparam/physics/radiative.F90 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/physics/radiative.F90 (revision 4197) +++ /dynamico_lmdz/simple_physics/phyparam/physics/radiative.F90 (revision 4198) @@ -2,112 +2,344 @@ IMPLICIT NONE SAVE + + PRIVATE + + PUBLIC :: mucorr, sw CONTAINS - - SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad) - + + PURE SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad) + + !======================================================================= + ! + ! Calcul of equivalent solar angle and and fraction of day whithout + ! diurnal cycle. + ! + ! parmeters : + ! ----------- + ! + ! Input : + ! ------- + ! npts number of points + ! pdeclin solar declinaison + ! plat(npts) latitude + ! phaut hauteur typique de l'atmosphere + ! prad rayon planetaire + ! + ! Output : + ! -------- + ! pmu(npts) equivalent cosinus of the solar angle + ! pfract(npts) fractionnal day + ! + !======================================================================= + + !----------------------------------------------------------------------- + ! + ! 0. Declarations : + ! ----------------- + + ! Arguments : + ! ----------- + INTEGER, INTENT(IN) :: npts + REAL, INTENT(IN) :: phaut, prad, pdeclin, plat(npts) + REAL, INTENT(OUT) :: pmu(npts), pfract(npts) + ! + ! Local variables : + ! ----------------- + INTEGER j + REAL z,cz,sz,tz,phi,cphi,sphi,tphi + REAL ap,a,t,b + REAL alph + + REAL, PARAMETER :: pi=2.*asin(1.) + + !----------------------------------------------------------------------- + +! print*,'npts,pdeclin',npts,pdeclin + z = pdeclin + cz = cos (z) + sz = sin (z) +! print*,'cz,sz',cz,sz + + DO j = 1, npts + + phi = plat(j) + cphi = cos(phi) + if (cphi.le.1.e-9) cphi=1.e-9 + sphi = sin(phi) + tphi = sphi / cphi + b = cphi * cz + t = -tphi * sz / cz + a = 1.0 - t*t + ap = a + + IF(t.eq.0.) then + t=0.5*pi + ELSE + IF (a.lt.0.) a = 0. + t = sqrt(a) / t + IF (t.lt.0.) then + t = -atan (-t) + pi + ELSE + t = atan(t) + ENDIF + ENDIF + + pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi + pfract(j) = t / pi + IF (ap .lt.0.) then + pmu(j) = sphi * sz + pfract(j) = 1.0 + ENDIF + IF (pmu(j).le.0.0) pmu(j) = 0.0 + pmu(j) = pmu(j) / pfract(j) + IF (pmu(j).eq.0.) pfract(j) = 0. + + END DO + + !----------------------------------------------------------------------- + ! correction de rotondite: + ! ------------------------ + + alph=phaut/prad + DO j=1,npts + ! !!!!!! + pmu(j)=sqrt(1224.*pmu(j)*pmu(j)+1.)/35. + ! $ (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j)) + END DO + + END SUBROUTINE mucorr + + PURE SUBROUTINE monGATHER(n,a,b,index) + INTEGER, INTENT(IN) :: n,index(n) + REAL, INTENT(IN) :: b(n) + REAL, INTENT(OUT) :: a(n) + INTEGER :: i + + DO i=1,n + a(i)=b(index(i)) + END DO + END SUBROUTINE monGATHER + + SUBROUTINE sw(ngrid,nlayer,ldiurn, & + coefvis,albedo, & + plevel,ps_rad,pmu,pfract,psolarf0, & + fsrfvis,dtsw, & + lwrite) + USE phys_const !======================================================================= ! -! Calcul of equivalent solar angle and and fraction of day whithout -! diurnal cycle. -! -! parmeters : -! ----------- -! -! Input : -! ------- -! npts number of points -! pdeclin solar declinaison -! plat(npts) latitude -! phaut hauteur typique de l'atmosphere -! prad rayon planetaire -! -! Output : -! -------- -! pmu(npts) equivalent cosinus of the solar angle -! pfract(npts) fractionnal day +! Rayonnement solaire en atmosphere non diffusante avec un +! coefficient d'absoprption gris. ! !======================================================================= - -!----------------------------------------------------------------------- -! -! 0. Declarations : -! ----------------- - -! Arguments : -! ----------- - INTEGER npts - REAL plat(npts),pmu(npts), pfract(npts) - REAL phaut,prad,pdeclin -! -! Local variables : -! ----------------- - INTEGER j - REAL pi - REAL z,cz,sz,tz,phi,cphi,sphi,tphi - REAL ap,a,t,b - REAL alph - -!----------------------------------------------------------------------- - - print*,'npts,pdeclin' - print*,npts,pdeclin - pi = 4. * atan(1.0) - print*,'PI=',pi - pi=2.*asin(1.) - print*,'PI=B',pi - z = pdeclin - cz = cos (z) - sz = sin (z) - print*,'cz,sz',cz,sz - - DO j = 1, npts - - phi = plat(j) - cphi = cos(phi) - if (cphi.le.1.e-9) cphi=1.e-9 - sphi = sin(phi) - tphi = sphi / cphi - b = cphi * cz - t = -tphi * sz / cz - a = 1.0 - t*t - ap = a - - IF(t.eq.0.) then - t=0.5*pi - ELSE - IF (a.lt.0.) a = 0. - t = sqrt(a) / t - IF (t.lt.0.) then - t = -atan (-t) + pi - ELSE - t = atan(t) +! +! declarations: +! ------------- +! +! +! arguments: +! ---------- +! + INTEGER ngrid,nlayer + REAL albedo(ngrid),coefvis + REAL pmu(ngrid),pfract(ngrid) + REAL plevel(ngrid,nlayer+1),ps_rad + REAL psolarf0 + REAL fsrfvis(ngrid),dtsw(ngrid,nlayer) + LOGICAL lwrite,ldiurn +! +! variables locales: +! ------------------ +! + + REAL zalb(ngrid),zmu(ngrid),zfract(ngrid) + REAL zplev(ngrid,nlayer+1) + REAL zflux(ngrid),zdtsw(ngrid,nlayer) + + INTEGER ig,l,nlevel,index(ngrid),ncount,igout + REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1) + REAL zfsrfref(ngrid) + REAL z1(ngrid) + REAL zu(ngrid,nlayer+1) + REAL tau0 + +!----------------------------------------------------------------------- +! 1. initialisations: +! ------------------- + + + nlevel=nlayer+1 + +!----------------------------------------------------------------------- +! Definitions des tableaux locaux pour les points ensoleilles: +! ------------------------------------------------------------ + + IF (ldiurn) THEN + ncount=0 + DO ig=1,ngrid + index(ig)=0 + ENDDO + DO ig=1,ngrid + IF(pfract(ig).GT.1.e-6) THEN + ncount=ncount+1 + index(ncount)=ig ENDIF - ENDIF - - pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi - pfract(j) = t / pi - IF (ap .lt.0.) then - pmu(j) = sphi * sz - pfract(j) = 1.0 - ENDIF - IF (pmu(j).le.0.0) pmu(j) = 0.0 - pmu(j) = pmu(j) / pfract(j) - IF (pmu(j).eq.0.) pfract(j) = 0. - - END DO - -!----------------------------------------------------------------------- -! correction de rotondite: + ENDDO + CALL monGATHER(ncount,zfract,pfract,index) + CALL monGATHER(ncount,zmu,pmu,index) + CALL monGATHER(ncount,zalb,albedo,index) + DO l=1,nlevel + CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index) + ENDDO + ELSE + ncount=ngrid + zfract(:)=pfract(:) + zmu(:)=pmu(:) + zalb(:)=albedo(:) + zplev(:,:)=plevel(:,:) + ENDIF + +!----------------------------------------------------------------------- +! calcul des profondeurs optiques integres depuis p=0: +! ---------------------------------------------------- + + tau0=-.5*log(coefvis) + +! calcul de la partie homogene de l'opacite + tau0=tau0/ps_rad + DO l=1,nlayer+1 + DO ig=1,ncount + zu(ig,l)=tau0*zplev(ig,l) + ENDDO + ENDDO + +!----------------------------------------------------------------------- +! 2. calcul de la transmission depuis le sommet de l'atmosphere: +! ----------------------------------------------------------- + + DO l=1,nlevel + DO ig=1,ncount + ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig)) + ENDDO + ENDDO + + IF (lwrite) THEN + igout=ncount/2+1 + PRINT* + PRINT*,'Diagnostique des transmission dans le spectre solaire' + PRINT*,'zfract, zmu, zalb' + PRINT*,zfract(igout), zmu(igout), zalb(igout) + PRINT*,'Pression, quantite d abs, transmission' + DO l=1,nlayer+1 + PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l) + ENDDO + ENDIF + +!----------------------------------------------------------------------- +! 3. taux de chauffage, ray. solaire direct: +! ------------------------------------------ + + DO l=1,nlayer + DO ig=1,ncount + zdtsw(ig,l)=g*psolarf0*zfract(ig)*zmu(ig)* & + (ztrdir(ig,l+1)-ztrdir(ig,l))/ & + (cpp*(zplev(ig,l)-zplev(ig,l+1))) + ENDDO + ENDDO + IF (lwrite) THEN + PRINT* + PRINT*,'Diagnostique des taux de chauffage solaires:' + PRINT*,' 1 taux de chauffage lie au ray. solaire direct' + DO l=1,nlayer + PRINT*,zdtsw(igout,l) + ENDDO + ENDIF + + +!----------------------------------------------------------------------- +! 4. calcul du flux solaire arrivant sur le sol: +! ---------------------------------------------- + + DO ig=1,ncount + z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1) + zflux(ig)=(1.-zalb(ig))*z1(ig) + zfsrfref(ig)= zalb(ig)*z1(ig) + ENDDO + IF (lwrite) THEN + PRINT* + PRINT*,'Diagnostique des taux de chauffage solaires:' + PRINT*,' 2 flux solaire net incident sur le sol' + PRINT*,zflux(igout) + ENDIF + + +!----------------------------------------------------------------------- +! 5.calcul des traansmissions depuis le sol, cas diffus: +! ------------------------------------------------------ + + DO l=1,nlevel + DO ig=1,ncount + ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66) + ENDDO + ENDDO + + IF (lwrite) THEN + PRINT* + PRINT*,'Diagnostique des taux de chauffage solaires' + PRINT*,' 3 transmission avec les sol' + PRINT*,'niveau transmission' + DO l=1,nlevel + PRINT*,l,ztrref(igout,l) + ENDDO + ENDIF + +!----------------------------------------------------------------------- +! 6.ajout a l'echauffement de la contribution du ray. sol. reflechit: +! ------------------------------------------------------------------- + + DO l=1,nlayer + DO ig=1,ncount + zdtsw(ig,l)=zdtsw(ig,l)+ & + g*zfsrfref(ig)*(ztrref(ig,l+1)-ztrref(ig,l))/ & + (cpp*(zplev(ig,l+1)-zplev(ig,l))) + ENDDO + ENDDO + +!----------------------------------------------------------------------- +! 10. sorites eventuelles: ! ------------------------ - alph=phaut/prad - DO j=1,npts -! !!!!!! - pmu(j)=sqrt(1224.*pmu(j)*pmu(j)+1.)/35. -! $ (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j)) - END DO - - END SUBROUTINE mucorr - - END MODULE radiative + IF (lwrite) THEN + PRINT* + PRINT*,'Diagnostique des taux de chauffage solaires:' + PRINT*,' 3 taux de chauffage total' + DO l=1,nlayer + PRINT*,zdtsw(igout,l) + ENDDO + ENDIF + + IF (ldiurn) THEN + fsrfvis(:)=0. + CALL monscatter(ncount,fsrfvis,index,zflux) + dtsw(:,:)=0. + DO l=1,nlayer + CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l)) + ENDDO + ELSE + print*,'NOT DIURNE' + fsrfvis(:)=zflux(:) + dtsw(:,:)=zdtsw(:,:) + ENDIF +! call dump2d(iim,jjm-1,zflux(2),'ZFLUX ') +! call dump2d(iim,jjm-1,fsrfvis(2),'FSRVIS ') +! call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') +! call dump2d(iim,jjm-1,pmu(2),'pmu ') +! call dump2d(iim,jjm-1,pfract(2),'pfract ') +! call dump2d(iim,jjm-1,albedo(2),'albedo ') +! call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') + + + END SUBROUTINE sw + +END MODULE radiative