source: dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90 @ 4233

Last change on this file since 4233 was 4229, checked in by dubos, 6 years ago

simple_physics : beautify code
File size: 16.0 KB
Line 
1MODULE phyparam_mod
2#include "use_logging.h"
3  USE callkeys
4  USE comgeomfi
5  IMPLICIT NONE
6  PRIVATE
7  SAVE
8
9  REAL :: nan ! initialized to NaN with adequate compiler options
10
11  REAL, PARAMETER :: pi=2*ASIN(1.), solarcst=1370., stephan=5.67e-08, &
12       ps_rad=1.e5, height_scale=10000., ref_temp=285., &
13       capcal_nosoil=1e5, tsoil_init=300.
14
15  ! internal state, written to / read from disk at checkpoint / restart
16  REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:), capcal(:), fluxgrd(:)
17  !$OMP THREADPRIVATE(tsurf, tsoil, capcal, fluxgrd)
18
19  ! precomputed arrays
20  REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:)
21  !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)
22
23  INTEGER :: icount
24  REAL    :: zday_last
25  !$OMP THREADPRIVATE( icount,zday_last)
26
27  PUBLIC :: phyparam
28
29CONTAINS
30
31  SUBROUTINE phyparam(ngrid,nlayer,             &
32       &            firstcall,lastcall,         &
33       &            rjourvrai,gmtime,ptimestep, &
34       &            pplev,pplay,pphi,           &
35       &            pu,pv,pt,                   &
36       &            pdu,pdv,pdt,pdpsrf)
37    USE phys_const, ONLY : g, rcp, r, unjours
38    USE surface,    ONLY : soil
39    USE turbulence, ONLY : vdif
40    USE convection, ONLY : convadj
41    USE writefield_mod, ONLY : writefield
42
43    !=======================================================================
44    !   Top routine of the physical parametrisations of the LMD
45    !   20 parameters GCM for planetary atmospheres.
46    !   It includes:
47    !   raditive transfer (long and shortwave) for CO2 and dust.
48    !   vertical turbulent mixing
49    !   convective adjsutment
50    !
51    !   author: Frederic Hourdin 15 / 10 /93
52    !=======================================================================
53
54    INTEGER, INTENT(IN) ::      &
55         ngrid,                 & ! Size of the horizontal grid.
56         nlayer                   ! Number of vertical layers.
57    LOGICAL, INTENT(IN) ::      &
58         firstcall,             & ! True at the first call
59         lastcall                 ! True at the last call
60    REAL, INTENT(IN)    ::      &
61         rjourvrai,             & ! Number of days counted from the North. Spring equinox
62         gmtime,                & ! time of the day in seconds
63         ptimestep,             & ! timestep (s)
64         pplev(ngrid,nlayer+1), & ! Pressure at interfaces between layers (pa)
65         pplay(ngrid,nlayer),   & ! Pressure at the middle of the layers (Pa)
66         pphi(ngrid,nlayer),    & ! Geopotential at the middle of the layers (m2s-2)
67         pu(ngrid,nlayer),      & ! u component of the wind (ms-1)
68         pv(ngrid,nlayer),      & ! v component of the wind (ms-1)
69         pt(ngrid,nlayer)         ! Temperature (K)
70    REAL, INTENT(OUT)   ::      & ! output : physical tendencies
71         pdu(ngrid,nlayer),     &
72         pdv(ngrid,nlayer),     &
73         pdt(ngrid,nlayer),     &
74         pdpsrf(ngrid)
75
76    !    Local variables :
77    REAL, DIMENSION(ngrid) :: mu0
78    INTEGER :: j,l,ig,nlevel,igout
79    !
80    REAL :: zday, zdtime
81    REAL zh(ngrid,nlayer),z1,z2
82    REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
83    REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)
84    REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
85    REAL zflubid(ngrid),zpmer(ngrid)
86    REAL zpopsk(ngrid,nlayer)
87    REAL zdum1(ngrid,nlayer)
88    REAL zdum2(ngrid,nlayer)
89    REAL zdum3(ngrid,nlayer)
90    REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
91    REAL zfluxsw(ngrid),zfluxlw(ngrid), fluxrad(ngrid)
92
93    WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid)
94    WRITELOG(*,*) 'nlayer',nlayer
95    LOG_DBG('phyparam')
96
97    IF (ngrid.NE.ngridmax) THEN
98       PRINT*,'STOP in inifis'
99       PRINT*,'Probleme de dimenesions :'
100       PRINT*,'ngrid     = ',ngrid
101       PRINT*,'ngridmax   = ',ngridmax
102       STOP
103    ENDIF
104
105    nlevel=nlayer+1
106    igout=ngrid/2+1
107    zday=rjourvrai+gmtime
108
109    !-----------------------------------------------------------------------
110    !    0. Allocate and initialize at first call
111    !    --------------------
112
113    IF(firstcall) THEN
114       !         zday_last=rjourvrai
115       zday_last=zday-ptimestep/unjours
116       CALL alloc(ngrid, nlayer)
117       CALL precompute
118       CALL coldstart(ngrid, ptimestep)
119    ENDIF
120
121    !-----------------------------------------------------------------------
122    !    1. Initialisations :
123    !    --------------------
124
125    icount=icount+1
126
127    pdv(:,:)  = 0.
128    pdu(:,:)  = 0.
129    pdt(:,:)  = 0.
130    pdpsrf(:) = 0.
131    zflubid(:)= 0.
132    zdtsrf(:) = 0.
133
134    !-----------------------------------------------------------------------
135    !   calcul du geopotentiel aux niveaux intercouches
136    !   ponderation des altitudes au niveau des couches en dp/p
137
138    DO l=1,nlayer
139       DO ig=1,ngrid
140          zzlay(ig,l)=pphi(ig,l)/g
141       ENDDO
142    ENDDO
143    DO ig=1,ngrid
144       zzlev(ig,1)=0.
145    ENDDO
146    DO l=2,nlayer
147       DO ig=1,ngrid
148          z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
149          z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
150          zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
151       ENDDO
152    ENDDO
153
154    !-----------------------------------------------------------------------
155    !   Transformation de la temperature en temperature potentielle
156    DO l=1,nlayer
157       DO ig=1,ngrid
158          zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp ! surface pressure is used as reference pressure
159          zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
160       ENDDO
161    ENDDO
162
163    !-----------------------------------------------------------------------
164    !    2. Compute radiative tendencies :
165    !    ---------------------------------------
166
167    IF(callrad) CALL radiative_tendencies(ngrid, igout, nlayer, &
168         gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
169         pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0)
170
171    !-----------------------------------------------------------------------
172    !    3. Vertical diffusion (turbulent mixing):
173    !    -----------------------------------------
174    !
175    IF(calldifv) THEN
176
177       DO ig=1,ngrid
178          zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
179       ENDDO
180
181       zdum1(:,:)=0.
182       zdum2(:,:)=0.
183
184       do l=1,nlayer
185          do ig=1,ngrid
186             zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)
187          enddo
188       enddo
189
190       CALL vdif(ngrid,nlayer,zday,        &
191            &        ptimestep,capcal,z0,       &
192            &        pplay,pplev,zzlay,zzlev,   &
193            &        pu,pv,zh,tsurf,emissiv,    &
194            &        zdum1,zdum2,zdum3,zflubid, &
195            &        zdufr,zdvfr,zdhfr,zdtsrfr, &
196            &        lverbose)
197
198       DO l=1,nlayer
199          DO ig=1,ngrid
200             pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
201             pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
202             pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
203          ENDDO
204       ENDDO
205
206       DO ig=1,ngrid
207          zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
208       ENDDO
209
210    ELSE
211       DO ig=1,ngrid
212          zdtsrf(ig)=zdtsrf(ig)+ &
213               &      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
214       ENDDO
215    ENDIF
216    !
217    !-----------------------------------------------------------------------
218    !   4. Dry convective adjustment:
219    !   -----------------------------
220
221    IF(calladj) THEN
222
223       DO l=1,nlayer
224          DO ig=1,ngrid
225             zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)
226          ENDDO
227       ENDDO
228
229       zdufr(:,:)=0.
230       zdvfr(:,:)=0.
231       zdhfr(:,:)=0.
232
233       CALL convadj(ngrid,nlayer,ptimestep, &
234            &                pplay,pplev,zpopsk, &
235            &                pu,pv,zh,           &
236            &                pdu,pdv,zdum1,      &
237            &                zdufr,zdvfr,zdhfr)
238
239       DO l=1,nlayer
240          DO ig=1,ngrid
241             pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
242             pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
243             pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
244          ENDDO
245       ENDDO
246
247    ENDIF
248
249    !-----------------------------------------------------------------------
250    !   On ajoute les tendances physiques a la temperature du sol:
251    !   ---------------------------------------------------------------
252
253    DO ig=1,ngrid
254       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
255    ENDDO
256
257    WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
258
259    !-----------------------------------------------------------------------
260    !   soil temperatures:
261    !   --------------------
262
263    IF (callsoil) THEN
264       CALL soil(ngrid,nsoilmx,.false.,inertie, &
265            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
266       IF(lverbose) THEN
267          WRITELOG(*,*) 'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
268          WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout), &
269               &        zdtsrf(igout), ptimestep
270          LOG_DBG('phyparam')
271       ENDIF
272    ENDIF
273
274    !-----------------------------------------------------------------------
275    !   sorties:
276    !   --------
277
278    WRITELOG(*,*) 'zday, zday_last ',zday,zday_last,icount
279    LOG_DBG('phyparam')
280
281    if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then
282       WRITELOG(*,*) 'zday, zday_last SORTIE ', zday, zday_last
283       LOG_INFO('phyparam')
284
285       zday_last=zday
286       !  Ecriture/extension de la coordonnee temps
287
288       do ig=1,ngridmax
289          zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*ref_temp))
290       enddo
291
292       call writefield('u','Vent zonal moy','m/s',pu)
293       call writefield('v','Vent meridien moy','m/s',pv)
294       call writefield('temp','Temperature','K',pt)
295       call writefield('geop','Geopotential','m2/s2',pphi)
296       call writefield('plev','plev','Pa',pplev(:,1:nlayer))
297
298       call writefield('du','du',' ',pdu)
299       call writefield('dv','du',' ',pdv)
300       call writefield('dt','du',' ',pdt)
301       call writefield('dtsw','dtsw',' ',zdtsw)
302       call writefield('dtlw','dtlw',' ',zdtlw)
303
304       call writefield('ts','Surface temper','K',tsurf)
305       call writefield('coslon','coslon',' ',coslon)
306       call writefield('sinlon','sinlon',' ',sinlon)
307       call writefield('coslat','coslat',' ',coslat)
308       call writefield('sinlat','sinlat',' ',sinlat)
309       call writefield('mu0','mu0',' ',mu0)
310       call writefield('alb','alb',' ',albedo)
311       call writefield('ps','Surface pressure','Pa',pplev(:,1))
312       call writefield('slp','Sea level pressure','Pa',zpmer)
313       call writefield('swsurf','SW surf','Pa',zfluxsw)
314       call writefield('lwsurf','LW surf','Pa',zfluxlw)
315
316    endif
317
318  END SUBROUTINE phyparam
319
320  SUBROUTINE radiative_tendencies(ngrid, igout, nlayer, &
321       gmtime, zdtime, zday, pplev, pplay, pt, &
322       pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0)
323    USE planet
324    USE phys_const,   ONLY : planet_rad, unjours
325    USE astronomy,    ONLY : orbite, solarlong
326    USE solar,        ONLY : solang, zenang, mucorr
327    USE radiative_sw, ONLY : sw
328    USE radiative_lw, ONLY : lw
329
330    INTEGER, INTENT(IN) :: ngrid, igout, nlayer
331    REAL, INTENT(IN)    :: gmtime, zdtime, zday, &
332         &                 pplev(ngrid,nlayer+1), pplay(ngrid, nlayer), &
333         &                 pt(ngrid, nlayer+1)
334    REAL, INTENT(INOUT) :: pdt(ngrid,nlayer)
335    REAL, INTENT(OUT)   :: zdtlw(ngrid,nlayer), & ! long-wave temperature tendency
336         &                 zfluxlw(ngrid),      & ! long-wave flux at surface
337         &                 zdtsw(ngrid,nlayer), & ! short-wave temperature tendency
338         &                 zfluxsw(ngrid),      & ! short-wave flux at surface
339         &                 fluxrad(ngrid),      & ! surface flux
340         mu0(ngrid)             ! ??
341    ! local variables
342    REAL :: fract(ngrid), dtrad(ngrid, nlayer)
343    REAL :: zls, zinsol, tsurf2, ztim1,ztim2,ztim3, dist_sol, declin
344    REAL :: zplanck(ngrid)
345    INTEGER :: ig, l
346
347    !    2.1 Insolation
348    !    --------------------------------------------------
349
350    CALL solarlong(zday,zls)
351    CALL orbite(zls,dist_sol,declin)
352
353    IF(diurnal) THEN
354       IF ( .TRUE. ) then
355          ztim1=SIN(declin)
356          ztim2=COS(declin)*COS(2.*pi*(zday-.5))
357          ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
358          CALL solang(ngrid,sinlon,coslon,sinlat,coslat, &
359               &         ztim1,ztim2,ztim3,                   &
360               &         mu0,fract)
361       ELSE
362          CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract)
363       ENDIF
364
365       IF(lverbose) THEN
366          WRITELOG(*,*) 'day, declin, sinlon,coslon,sinlat,coslat'
367          WRITELOG(*,*) zday, declin,                       &
368               &        sinlon(igout),coslon(igout),  &
369               &        sinlat(igout),coslat(igout)
370          LOG_DBG('radiative_tendencies')
371       ENDIF
372    ELSE
373       WRITELOG(*,*) 'declin,ngrid,planet_rad', declin, ngrid, planet_rad
374       LOG_DBG('radiative_tendencies')
375       CALL mucorr(ngrid,declin,lati,mu0,fract,height_scale,planet_rad)
376    ENDIF
377
378    zinsol=solarcst/(dist_sol*dist_sol)
379
380    !    2.2 Radiative tendencies and fluxes:
381    !    --------------------------------------------------
382
383    CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, &
384         &              pplev,ps_rad,                 &
385         &              mu0,fract,zinsol,             &
386         &              zfluxsw,zdtsw,                &
387         &              lverbose)
388
389    CALL lw(ngrid,nlayer,coefir,emissiv, &
390         &             pplev,ps_rad,tsurf,pt, &
391         &             zfluxlw,zdtlw,         &
392         &             lverbose)
393
394    !    2.4 surface fluxes
395    !    ------------------------------
396
397    DO ig=1,ngrid
398       fluxrad(ig)=emissiv(ig)*zfluxlw(ig) &
399            &         +zfluxsw(ig)*(1.-albedo(ig))
400       tsurf2 = tsurf(ig)*tsurf(ig)
401       zplanck(ig)=emissiv(ig)*stephan*tsurf2*tsurf2
402       fluxrad(ig)=fluxrad(ig)-zplanck(ig)
403    ENDDO
404
405    !    2.5 Temperature tendencies
406    !    --------------------------
407
408    DO l=1,nlayer
409       DO ig=1,ngrid
410          dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
411          pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
412       ENDDO
413    ENDDO
414
415    IF(lverbose) THEN
416       WRITELOG(*,*) 'Diagnostics for radiation'
417       WRITELOG(*,*) 'albedo, emissiv, mu0,fract,Frad,Planck'
418       WRITELOG(*,*) albedo(igout),emissiv(igout),mu0(igout), &
419            &        fract(igout),                       &
420            &        fluxrad(igout),zplanck(igout)
421       WRITELOG(*,*) 'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
422       WRITELOG(*,*) 'unjours',unjours
423       DO l=1,nlayer
424          WRITELOG(*,'(3f15.5,2e15.2)') pt(igout,l),     &
425               &         pplay(igout,l),pplev(igout,l),  &
426               &         zdtsw(igout,l),zdtlw(igout,l)
427       ENDDO
428       LOG_DBG('radiative_tendencies')
429    ENDIF
430
431  END SUBROUTINE radiative_tendencies
432
433  SUBROUTINE alloc(ngrid, nlayer)
434    USE astronomy, ONLY : iniorbit
435    USE surface, ONLY : zc,zd
436    INTEGER, INTENT(IN) :: ngrid, nlayer
437    LOGICAL, PARAMETER :: firstcall=.TRUE.
438    ! allocate arrays for internal state
439    ALLOCATE(tsurf(ngrid))
440    ALLOCATE(tsoil(ngrid,nsoilmx))
441    ! we could avoid the arrays below with a different implementation of surface / radiation / turbulence coupling
442    ALLOCATE(capcal(ngrid),fluxgrd(ngrid))
443    ALLOCATE(zc(ngrid,nsoilmx),zd(ngrid,nsoilmx))
444    ! allocate precomputed arrays
445    ALLOCATE(rnatur(ngrid))
446    ALLOCATE(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
447    CALL iniorbit
448  END SUBROUTINE alloc
449
450  SUBROUTINE precompute
451    USE surface
452    ! precompute time-independent arrays
453    rnatur(:)  = 1.
454    emissiv(:) = (1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter
455    inertie(:) = (1.-rnatur(:))*I_mer+rnatur(:)*I_ter
456    albedo(:)  = (1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter
457    z0(:)      = (1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter
458  END SUBROUTINE precompute
459
460  SUBROUTINE coldstart(ngrid, ptimestep)
461    ! create internal state to start a run without a restart file
462    USE surface, ONLY : soil
463    INTEGER, INTENT(IN) :: ngrid
464    REAL, INTENT(iN)    :: ptimestep
465    tsurf(:)   = tsoil_init
466    tsoil(:,:) = tsoil_init
467    icount=0
468    IF(callsoil) THEN
469       ! initializes zc, zd, capcal, fluxgrd
470       CALL soil(ngrid,nsoilmx,.TRUE.,inertie, &
471            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
472    ELSE
473       WRITELOG(*,*) 'WARNING!!! Thermal conduction in the soil turned off'
474       LOG_WARN('coldstart')
475       capcal(:)  = capcal_nosoil
476       fluxgrd(:) = 0.
477    ENDIF
478  END SUBROUTINE coldstart
479
480END MODULE phyparam_mod
Note: See TracBrowser for help on using the repository browser.