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

Last change on this file since 4240 was 4236, checked in by dubos, 5 years ago

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