Changeset 4198 for dynamico_lmdz

Timestamp:

Dec 20, 2019, 12:00:02 PM (5 years ago)

Author:

dubos

Message:

simple_physics : cleanup radiative transfer

Location:

dynamico_lmdz/simple_physics/phyparam

Files:

• param/gather.F (deleted)
• param/phyparam.F (modified) (1 diff)
• param/sw.F (deleted)
• physics/radiative.F90 (modified) (1 diff)

dynamico_lmdz/simple_physics/phyparam/param/phyparam.F

@@ -14 +14 @@
       USE astronomy
       USE vdif_mod, ONLY : vdif
-      USE radiative, ONLY : mucorr
+      USE radiative, ONLY : mucorr, sw
 c     
       IMPLICIT NONE

dynamico_lmdz/simple_physics/phyparam/physics/radiative.F90

@@ -2 +2 @@
   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