source: dynamico_lmdz/simple_physics/phyparam/physics/radiative_lw.F90 @ 4221

Last change on this file since 4221 was 4199, checked in by dubos, 5 years ago

simple_physics : cleanup radiative transfer
File size: 6.9 KB
Line 
1MODULE radiative_lw
2
3#include "use_logging.h"
4
5  IMPLICIT NONE
6  SAVE
7 
8  PRIVATE
9 
10  PUBLIC :: lw
11 
12  LOGICAL, PARAMETER :: lstrong=.TRUE.
13  REAL, PARAMETER :: stephan=5.67e-08
14 
15CONTAINS
16 
17  SUBROUTINE lw(ngrid,nlayer,coefir,emissiv, &
18       pp,ps_rad,ptsurf,pt,              &
19       pfluxir,pdtlw,                    &
20       lwrite)
21    USE phys_const, ONLY : cpp, g
22    !=======================================================================
23    !
24    !   calcul de l'evolution de la temperature sous l'effet du rayonnement
25    !   infra-rouge.
26    !   Pour simplifier, les transmissions sont precalculees et ne
27    !   dependent que de l'altitude.
28    !
29    !   arguments:
30    !   ----------
31    !
32    !   entree:
33    !   -------
34    !      ngrid             nombres de points de la grille horizontale
35    !      nlayer              nombre de couches
36    !      ptsurf(ngrid)     temperature de la surface
37    !      pt(ngrid,nlayer)    temperature des couches
38    !      pp(ngrid,nlayer+1)  pression entre les couches
39    !      lwrite            variable logique pour sorties
40    !
41    !   sortie:
42    !   -------
43    !      pdtlw(ngrid,nlayer) taux de refroidissement
44    !      pfluxir(ngrid)    flux infrarouge sur le sol
45    !
46    !=======================================================================
47   
48    !   declarations:
49    !   -------------
50   
51    !   arguments:
52    !   ----------
53   
54    INTEGER, INTENT(IN)  ::  ngrid,nlayer
55    REAL,    INTENT(IN)  :: coefir,emissiv(ngrid),ps_rad
56    REAL,    INTENT(IN)  ::  ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)
57    REAL,    INTENT(OUT) ::  pdtlw(ngrid,nlayer),pfluxir(ngrid)
58    LOGICAL, INTENT(IN)  :: lwrite
59   
60    !   variables locales:
61    !   ------------------
62   
63    INTEGER nlevel,ilev,ig,i,il
64    REAL zplanck(ngrid,nlayer+1),zcoef
65    REAL zfluxup(ngrid,nlayer+1),zfluxdn(ngrid,nlayer+1)
66    REAL zflux(ngrid,nlayer+1)
67    REAL zlwtr1(ngrid),zlwtr2(ngrid)
68    REAL zup(ngrid,nlayer+1),zdup(ngrid)
69
70    CHARACTER(:), PARAMETER :: tag='rad/lw'
71    !-----------------------------------------------------------------------
72    !   initialisations:
73    !   ----------------
74   
75    nlevel=nlayer+1
76   
77    !-----------------------------------------------------------------------
78    !   2. calcul des quantites d'absorbants:
79    !   -------------------------------------
80   
81    !   absorption forte
82    IF(lstrong) THEN
83       DO ilev=1,nlevel
84          DO ig=1,ngrid
85             zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g)
86          ENDDO
87       ENDDO
88       IF(lwrite) THEN
89          DO ilev=1,nlayer
90             WRITELOG(*,*) ' up(',ilev,')  =  ',zup(ngrid/2+1,ilev)
91          ENDDO
92          LOG_DBG(tag)
93       ENDIF
94       zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))
95       
96       !   absorption faible
97    ELSE
98       DO ilev=1,nlevel
99          DO ig=1,ngrid
100             zup(ig,ilev)=pp(ig,ilev)
101          ENDDO
102       ENDDO
103       zcoef=-log(coefir)/ps_rad
104    ENDIF
105   
106   
107    !-----------------------------------------------------------------------
108    !   2. calcul de la fonction de corps noir:
109    !   ---------------------------------------
110   
111    DO ilev=1,nlayer
112       DO ig=1,ngrid
113          zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)
114          zplanck(ig,ilev)=stephan* &
115               zplanck(ig,ilev)*zplanck(ig,ilev)
116       ENDDO
117    ENDDO
118   
119    !-----------------------------------------------------------------------
120    !   4. flux descendants:
121    !   --------------------
122   
123    DO ilev=1,nlayer
124       DO ig=1,ngrid
125          zfluxdn(ig,ilev)=0.
126       ENDDO
127       DO ig=1,ngrid
128          zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)
129       ENDDO
130       CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
131       
132       DO il=nlayer,ilev,-1
133          zlwtr2(:)=zlwtr1(:)
134          DO ig=1,ngrid
135             zdup(ig)=zup(ig,ilev)-zup(ig,il)
136          ENDDO
137          CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
138          DO ig=1,ngrid
139             zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+ &
140                  zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
141          ENDDO
142       ENDDO
143    ENDDO
144   
145    DO ig=1,ngrid
146       zfluxdn(ig,nlevel)=0.
147       pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
148    ENDDO
149   
150    DO ig=1,ngrid
151       zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
152       zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1) &
153            +(1.-emissiv(ig))*zfluxdn(ig,1)
154    ENDDO
155   
156    !-----------------------------------------------------------------------
157    !   3. flux montants:
158    !   ------------------
159   
160    DO ilev=1,nlayer
161       DO ig=1,ngrid
162          zdup(ig)=zup(ig,1)-zup(ig,ilev+1)
163       ENDDO
164       CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
165       DO ig=1,ngrid
166          zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)
167       ENDDO
168       DO il=1,ilev
169          zlwtr2(:)=zlwtr1(:)
170          DO ig=1,ngrid
171             zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)
172          ENDDO
173          CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
174          DO ig=1,ngrid
175             zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+ &
176                  zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
177          ENDDO
178       ENDDO
179       
180    ENDDO
181
182    !-----------------------------------------------------------------------
183    !   5. calcul des flux nets:
184    !   ------------------------
185   
186    DO ilev=1,nlevel
187       DO ig=1,ngrid
188          zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)
189       ENDDO
190    ENDDO
191   
192    !-----------------------------------------------------------------------
193    !   6. Calcul des taux de refroidissement:
194    !   --------------------------------------
195   
196    DO ilev=1,nlayer
197       DO ig=1,ngrid
198          pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))* &
199               g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))
200       ENDDO
201    ENDDO
202   
203    !-----------------------------------------------------------------------
204    !   10. sorties eventuelles:
205    !   ------------------------
206   
207    IF (lwrite) THEN       
208       WRITELOG(*,*) 'Diagnostique rayonnement thermique'
209       WRITELOG(*,*) 'temperature     ', &
210          'flux montant    flux desc.     taux de refroid.'
211       i=ngrid/2+1
212       WRITELOG(6,'(4e18.4)') ptsurf(i)
213       DO ilev=1,nlayer
214          WRITELOG(6,'(i4,4e18.4)') ilev,pt(i,ilev), &
215               zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
216       ENDDO
217       WRITELOG(6,'(4e18.4)') zfluxup(i,nlevel),zfluxdn(i,nlevel)       
218       LOG_DBG(tag)
219    ENDIF
220
221!-----------------------------------------------------------------------
222   
223  END SUBROUTINE lw
224
225  PURE SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)
226    INTEGER, INTENT(IN) :: ngrid
227    REAL,    INTENT(IN) :: coef
228    LOGICAL, INTENT(IN) :: lstrong
229    REAL,    INTENT(IN) :: dup(ngrid)
230    REAL,    INTENT(OUT) :: transm(ngrid)
231    INTEGER ig
232    IF(lstrong) THEN
233       DO ig=1,ngrid
234          transm(ig)=exp(-coef*sqrt(dup(ig)))
235       ENDDO
236    ELSE
237       DO ig=1,ngrid
238          transm(ig)=exp(-coef*dup(ig))
239       ENDDO
240    ENDIF
241   
242  END SUBROUTINE lwtr
243 
244END MODULE radiative_lw
Note: See TracBrowser for help on using the repository browser.