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

Last change on this file was 4248, checked in by dubos, 5 years ago

simple_physics : output theta
File size: 13.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, PARAMETER :: ref_temp=285., capcal_nosoil=1e5, tsoil_init=300.
10
11  INTEGER :: icount
12  REAL    :: zday_last
13  !$OMP THREADPRIVATE( icount,zday_last)
14
15  PUBLIC :: phyparam
16
17CONTAINS
18
19  SUBROUTINE phyparam(ngrid,nlayer,             &
20       &            firstcall,lastcall,         &
21       &            rjourvrai,gmtime,ptimestep, &
22       &            pplev,pplay,pphi,           &
23       &            pu,pv,pt,                   &
24       &            pdu,pdv,pdt,pdpsrf) BIND(C, name='phyparam_phyparam')
25    USE phys_const, ONLY : g, rcp, r, unjours
26    USE soil_mod, ONLY : soil_forward, soil_backward
27    USE soil_mod, ONLY : z0, inertie, emissiv, albedo  ! precomputed
28    USE soil_mod, ONLY : tsurf, tsoil ! state variables
29    USE turbulence, ONLY : vdif
30    USE convection, ONLY : convadj
31    USE radiative_mod,  ONLY : radiative_tendencies
32    USE writefield_mod, ONLY : writefield
33
34    !=======================================================================
35    !   Top routine of the physical parametrisations of the LMD
36    !   20 parameters GCM for planetary atmospheres.
37    !   It includes:
38    !   radiative transfer (long and shortwave) for CO2 and dust.
39    !   vertical turbulent mixing
40    !   convective adjsutment
41    !   heat diffusion in the soil
42    !
43    !   author: Frederic Hourdin 15 / 10 /93
44    !=======================================================================
45
46    INTEGER, INTENT(IN), VALUE :: &
47         ngrid,                 & ! Size of the horizontal grid.
48         nlayer                   ! Number of vertical layers.
49    LOGICAL, INTENT(IN), VALUE  :: &
50         firstcall,             & ! True at the first call
51         lastcall                 ! True at the last call
52    REAL, INTENT(IN), VALUE     ::      &
53         rjourvrai,             & ! Number of days counted from the North. Spring equinox
54         gmtime,                & ! time of the day in seconds
55         ptimestep                ! timestep (s)
56    REAL, INTENT(IN) :: &
57         pplev(ngrid,nlayer+1), & ! Pressure at interfaces between layers (pa)
58         pplay(ngrid,nlayer),   & ! Pressure at the middle of the layers (Pa)
59         pphi(ngrid,nlayer),    & ! Geopotential at the middle of the layers (m2s-2)
60         pu(ngrid,nlayer),      & ! u component of the wind (ms-1)
61         pv(ngrid,nlayer),      & ! v component of the wind (ms-1)
62         pt(ngrid,nlayer)         ! Temperature (K)
63    REAL, INTENT(OUT)   ::      & ! output : physical tendencies
64         pdu(ngrid,nlayer),     &
65         pdv(ngrid,nlayer),     &
66         pdt(ngrid,nlayer),     &
67         pdpsrf(ngrid)
68
69    !    Local variables :
70    REAL :: zh(ngrid,nlayer),      & ! potential temperature
71         &  zpopsk(ngrid,nlayer),  & ! Exner function
72         &  zzlev(ngrid,nlayer+1), & ! altitude of interfaces
73         &  zzlay(ngrid,nlayer),   & ! altitude of full levels
74         &  fluxrad(ngrid),        & ! radiative flux at surface
75         &  zc(ngrid, nsoilmx),    & ! LU coefficients for soil implicit solve
76         &  zd(ngrid, nsoilmx),    &
77         &  fluxgrd(ngrid),        & ! heat flux from deep soil
78         &  capcal(ngrid),         & ! effective heat capacity of soil
79         &  zdufr(ngrid,nlayer),   & ! partial tendencies for zonal velocity,
80         &  zdvfr(ngrid,nlayer),   & !   meridional velocity,
81         &  zdhfr(ngrid,nlayer),   & !   potential temperature,
82         &  zdtsrfr(ngrid),        & !   surface temperature
83         &  zdtsrf(ngrid),         & ! total tendency of surface temperature
84         &  zflubid(ngrid),        & ! radiative + deep soil fluxes
85         &  zpmer(ngrid)            ! sea-level pressure
86    REAL zdum1(ngrid,nlayer)
87    REAL zdum2(ngrid,nlayer)
88    REAL zdum3(ngrid,nlayer)
89
90    INTEGER :: j,l,ig,nlevel,igout
91    LOGICAL :: lwrite
92    REAL    :: zday, zdtime, z1, z2
93
94    WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid)
95    WRITELOG(*,*) 'nlayer',nlayer
96    LOG_DBG('phyparam')
97
98    IF (ngrid.NE.ngridmax) THEN
99       PRINT*,'STOP in inifis'
100       PRINT*,'Probleme de dimensions :'
101       PRINT*,'ngrid     = ',ngrid
102       PRINT*,'ngridmax   = ',ngridmax
103       STOP
104    ENDIF
105
106    nlevel=nlayer+1
107    igout=ngrid/2+1
108    zday=rjourvrai+gmtime
109
110    !-----------------------------------------------------------------------
111    !    0. Allocate and initialize at first call
112    !    --------------------
113
114    IF(firstcall) THEN
115       !         zday_last=rjourvrai
116       zday_last=zday-ptimestep/unjours
117       CALL alloc(ngrid, nlayer)
118       CALL precompute
119       CALL coldstart(ngrid)
120    ENDIF
121
122    !-----------------------------------------------------------------------
123    !    1. Initialisations :
124    !    --------------------
125
126    icount=icount+1
127
128    pdv(:,:)  = 0.
129    pdu(:,:)  = 0.
130    pdt(:,:)  = 0.
131    pdpsrf(:) = 0.
132    zflubid(:)= 0.
133    zdtsrf(:) = 0.
134
135    IF(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) THEN
136       WRITELOG(*,*) 'zday, zday_last SORTIE ', zday, zday_last
137       LOG_INFO('phyparam')
138       zday_last=zday
139       lwrite = .TRUE.
140    END IF
141
142    !-----------------------------------------------------------------------
143    !   calcul du geopotentiel aux niveaux intercouches
144    !   ponderation des altitudes au niveau des couches en dp/p
145
146    DO l=1,nlayer
147       DO ig=1,ngrid
148          zzlay(ig,l)=pphi(ig,l)/g
149       ENDDO
150    ENDDO
151    DO ig=1,ngrid
152       zzlev(ig,1)=0.
153    ENDDO
154    DO l=2,nlayer
155       DO ig=1,ngrid
156          z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
157          z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
158          zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
159       ENDDO
160    ENDDO
161
162    !-----------------------------------------------------------------------
163    !   Transformation de la temperature en temperature potentielle
164    DO l=1,nlayer
165       DO ig=1,ngrid
166          zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp ! surface pressure is used as reference pressure
167          zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
168       ENDDO
169    ENDDO
170
171    !-------------------------------------------------------------
172    !  soil temperatures : 1st half of implicit time integration
173    !  forward sweep from deep ground to surface
174    !  yields LU coefficients zc,zd and capcal, fluxgrd
175    !   ----------------------------------------------------------
176
177    IF (callsoil) THEN
178       CALL soil_forward(ngrid,nsoilmx, ptimestep, inertie, tsurf, tsoil, &
179            &            zc, zd, capcal, fluxgrd)
180
181       !       CALL soil_new(ngrid,nsoilmx,ptimestep,inertie, &
182       !            tsurf, tsoil, capcal,fluxgrd)
183       !       CALL soil(ngrid,nsoilmx,.false.,inertie, &
184       !            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
185    ELSE
186       capcal(:)  = capcal_nosoil
187       fluxgrd(:) = 0.
188    END IF
189
190    IF(lverbose) THEN
191       WRITELOG(*,*) 'Surface Heat capacity, conduction Flux, Ts'
192       WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout)
193       LOG_DBG('phyparam')
194    END IF
195
196    !-----------------------------------------------------------------------
197    !    2. Compute radiative tendencies :
198    !    ---------------------------------------
199
200    IF(callrad) CALL radiative_tendencies(lwrite, ngrid, igout, nlayer, &
201         gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
202         pdt, fluxrad)
203
204    !-----------------------------------------------------------------------
205    !    3. Vertical diffusion (turbulent mixing):
206    ! Kz is computed then vertical diffusion is integrated in time implicitly
207    ! using a linear relationship between surface heat flux and air temperature
208    ! in lowest level (Robin-type BC)
209    !    -------------------------------------------------------------------
210    !
211    IF(calldifv) THEN
212
213       DO ig=1,ngrid
214          zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
215       ENDDO
216
217       zdum1(:,:)=0.
218       zdum2(:,:)=0.
219
220       do l=1,nlayer
221          do ig=1,ngrid
222             zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)
223          enddo
224       enddo
225
226       CALL vdif(ngrid,nlayer,zday,        &
227            &        ptimestep,capcal,z0,       &
228            &        pplay,pplev,zzlay,zzlev,   &
229            &        pu,pv,zh,tsurf,emissiv,    &
230            &        zdum1,zdum2,zdum3,zflubid, &
231            &        zdufr,zdvfr,zdhfr,zdtsrfr, &
232            &        lverbose)
233
234       DO l=1,nlayer
235          DO ig=1,ngrid
236             pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
237             pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
238             pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
239          ENDDO
240       ENDDO
241
242       DO ig=1,ngrid
243          zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
244       ENDDO
245
246    ELSE
247       DO ig=1,ngrid
248          zdtsrf(ig)=zdtsrf(ig)+ &
249               &      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
250       ENDDO
251    ENDIF
252
253    !-------------------------------------------------------------
254    !   soil temperatures : 2nd half of implicit time integration
255    !   using updated tsurf as input
256    !   ----------------------------------------------------------
257
258    DO ig=1,ngrid
259       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
260    ENDDO
261
262    WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
263
264    IF (callsoil) THEN
265       CALL soil_backward(ngrid,nsoilmx, zc,zd, tsurf,tsoil)
266       IF(lverbose) THEN
267          WRITELOG(*,*) 'Surface Ts, dTs, dt'
268          WRITELOG(*,*) tsurf(igout), zdtsrf(igout), ptimestep
269          LOG_DBG('phyparam')
270       ENDIF
271    END IF
272
273
274    !
275    !-----------------------------------------------------------------------
276    !   4. Dry convective adjustment:
277    !   -----------------------------
278
279    IF(calladj) THEN
280
281       DO l=1,nlayer
282          DO ig=1,ngrid
283             zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)
284          ENDDO
285       ENDDO
286
287       zdufr(:,:)=0.
288       zdvfr(:,:)=0.
289       zdhfr(:,:)=0.
290
291       CALL convadj(ngrid,nlayer,ptimestep, &
292            &                pplay,pplev,zpopsk, &
293            &                pu,pv,zh,           &
294            &                pdu,pdv,zdum1,      &
295            &                zdufr,zdvfr,zdhfr)
296
297       DO l=1,nlayer
298          DO ig=1,ngrid
299             pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
300             pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
301             pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
302          ENDDO
303       ENDDO
304
305    ENDIF
306
307    !-----------------------------------------------------------------------
308    !   sorties:
309    !   --------
310
311    WRITELOG(*,*) 'zday, zday_last ',zday,zday_last,icount
312    LOG_DBG('phyparam')
313
314    IF(lwrite) THEN
315
316       do ig=1,ngridmax
317          zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*ref_temp))
318       enddo
319
320       call writefield('u','Vent zonal moy','m/s',pu)
321       call writefield('v','Vent meridien moy','m/s',pv)
322       call writefield('temp','Temperature','K',pt)
323       call writefield('theta','Potential temperature','K',zh)
324       call writefield('geop','Geopotential','m2/s2',pphi)
325       call writefield('plev','plev','Pa',pplev(:,1:nlayer))
326
327       call writefield('du','du',' ',pdu)
328       call writefield('dv','dv',' ',pdv)
329       call writefield('dt','dt',' ',pdt)
330
331       call writefield('ts','Surface temper','K',tsurf)
332       call writefield('coslon','coslon',' ',coslon)
333       call writefield('sinlon','sinlon',' ',sinlon)
334       call writefield('coslat','coslat',' ',coslat)
335       call writefield('sinlat','sinlat',' ',sinlat)
336       call writefield('alb','alb',' ',albedo)
337       call writefield('ps','Surface pressure','Pa',pplev(:,1))
338       call writefield('slp','Sea level pressure','Pa',zpmer)
339    END IF
340
341  END SUBROUTINE phyparam
342
343  SUBROUTINE alloc(ngrid, nlayer) BIND(C, name='phyparam_alloc')
344    !$cython header void phyparam_alloc(int, int);
345    !$cython wrapper def alloc(ngrid, nlayer) : phy.phyparam_alloc(ngrid, nlayer)
346    USE astronomy, ONLY : iniorbit
347    USE soil_mod
348    INTEGER, INTENT(IN), VALUE :: ngrid, nlayer
349    ! allocate precomputed arrays
350    ALLOCATE(rnatur(ngrid), albedo(ngrid), emissiv(ngrid))
351    ALLOCATE(z0(ngrid),inertie(ngrid))
352    ! allocate arrays for internal state
353    ALLOCATE(tsurf(ngrid))
354    ALLOCATE(tsoil(ngrid,nsoilmx))
355    CALL iniorbit
356  END SUBROUTINE alloc
357
358  SUBROUTINE precompute() BIND(C, name='phyparam_precompute')
359    !$cython header void phyparam_precompute();
360    !$cython wrapper def precompute() : phy.phyparam_precompute()
361    ! precompute time-independent arrays
362    USE soil_mod
363    rnatur(:)  = 1.
364    inertie(:) = (1.-rnatur(:))*I_mer+rnatur(:)*I_ter
365    z0(:)      = (1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter
366    emissiv(:) = (1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter
367    albedo(:)  = (1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter
368  END SUBROUTINE precompute
369
370  SUBROUTINE coldstart(ngrid) BIND(C, name='phyparam_coldstart')
371    !$cython header void phyparam_coldstart(int);
372    !$cython wrapper def coldstart (ngrid): phy.phyparam_coldstart(ngrid)
373    ! create internal state to start a run without a restart file
374    USE soil_mod, ONLY : tsurf, tsoil
375    INTEGER, INTENT(IN), VALUE :: ngrid
376    tsurf(:)   = tsoil_init
377    tsoil(:,:) = tsoil_init
378    icount=0
379    IF(.NOT. callsoil) THEN
380       WRITELOG(*,*) 'WARNING!!! Thermal conduction in the soil turned off'
381       LOG_WARN('coldstart')
382    ENDIF
383  END SUBROUTINE coldstart
384
385END MODULE phyparam_mod
Note: See TracBrowser for help on using the repository browser.