source: trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90 @ 797

Last change on this file since 797 was 787, checked in by aslmd, 13 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

  • nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this
File size: 10.2 KB
RevLine 
[787]1subroutine moistadj(ngrid, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
[135]2
[728]3  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, Psat_water, Lcpdqsat_water
[787]4  USE tracer_h
[135]5
6  implicit none
7
8
9!=====================================================================
10!     
11!     Purpose
12!     -------
13!     Calculates moist convective adjustment by the method of Manabe.
14!     
15!     Authors
16!     -------
17!     Adapted from the LMDTERRE code by R. Wordsworth (2010)
18!     Original author Z. X. Li (1993)
19!     
20!=====================================================================
21
22#include "dimensions.h"
23#include "dimphys.h"
24#include "comcstfi.h"
25
[787]26      INTEGER ngrid, nq
[135]27
[787]28      REAL pt(ngrid,nlayermx)            ! temperature (K)
29      REAL pq(ngrid,nlayermx,nq)       ! tracer (kg/kg)
30      REAL pdq(ngrid,nlayermx,nq)
[135]31
[787]32      REAL pdqmana(ngrid,nlayermx,nq)  ! tendency of tracers (kg/kg.s-1)
33      REAL pdtmana(ngrid,nlayermx)       ! temperature increment
34
[728]35!     local variables
[787]36      REAL zt(ngrid,nlayermx)      ! temperature (K)
37      REAL zq(ngrid,nlayermx)      ! humidite specifique (kg/kg)
38      REAL pplev(ngrid,nlayermx+1) ! pression a inter-couche (Pa)
39      REAL pplay(ngrid,nlayermx)   ! pression au milieu de couche (Pa)
[135]40
[787]41      REAL d_t(ngrid,nlayermx)     ! temperature increment
42      REAL d_q(ngrid,nlayermx)     ! incrementation pour vapeur d'eau
43      REAL d_ql(ngrid,nlayermx)    ! incrementation pour l'eau liquide
44      REAL rneb(ngrid,nlayermx) ! cloud fraction
[135]45      REAL ptimestep
46
47!      REAL t_coup
48!      PARAMETER (t_coup=234.0)
49      REAL seuil_vap
50      PARAMETER (seuil_vap=1.0E-10)
51
52!     Local variables
[728]53      INTEGER i, k, iq, nn
[773]54      INTEGER, PARAMETER :: niter = 1
[135]55      INTEGER k1, k1p, k2, k2p
[787]56      LOGICAL itest(ngrid)
57      REAL delta_q(ngrid, nlayermx)
[135]58      REAL cp_new_t(nlayermx)
59      REAL cp_delta_t(nlayermx)
60      REAL v_cptj(nlayermx), v_cptjk1, v_ssig
[787]61      REAL v_cptt(ngrid,nlayermx), v_p, v_t, zpsat
62      REAL zqs(ngrid,nlayermx), zdqs(ngrid,nlayermx)
63      REAL zq1(ngrid), zq2(ngrid)
64      REAL gamcpdz(ngrid,2:nlayermx)
[135]65      REAL zdp, zdpm
66
67      REAL zsat ! super-saturation
68      REAL zflo ! flotabilite
69
[787]70      REAL local_q(ngrid,nlayermx),local_t(ngrid,nlayermx)
[135]71
72      REAL zdelta, zcor, zcvm5
73
74      REAL dEtot, dqtot, masse ! conservation diagnostics
75      real dL1tot, dL2tot
76
77!     Indices of water vapour and water ice tracers
78      INTEGER,SAVE :: i_h2o=0  ! water vapour
79      INTEGER,SAVE :: i_ice=0  ! water ice
80
81      LOGICAL firstcall
82      SAVE firstcall
83
84      DATA firstcall /.TRUE./
85
86      IF (firstcall) THEN
87
88         i_h2o=igcm_h2o_vap
89         i_ice=igcm_h2o_ice
90       
91         write(*,*) "rain: i_ice=",i_ice
92         write(*,*) "      i_h2o=",i_h2o
93
94         firstcall = .FALSE.
95      ENDIF
96
97!     GCM -----> subroutine variables
[787]98      zq(1:ngrid,1:nlayermx)    = pq(1:ngrid,1:nlayermx,i_h2o)+ pdq(1:ngrid,1:nlayermx,i_h2o)*ptimestep
99      zt(1:ngrid,1:nlayermx)    = pt(1:ngrid,1:nlayermx)
100      pdqmana(1:ngrid,1:nlayermx,1:nq)=0.0
[728]101
[135]102      DO k = 1, nlayermx
[787]103       DO i = 1, ngrid
[728]104         if(zq(i,k).lt.0.)then
105            zq(i,k)=0.0
[135]106         endif
[728]107       ENDDO
[135]108      ENDDO
[728]109     
[787]110      local_q(1:ngrid,1:nlayermx) = zq(1:ngrid,1:nlayermx)
111      local_t(1:ngrid,1:nlayermx) = zt(1:ngrid,1:nlayermx)
112      rneb(1:ngrid,1:nlayermx) = 0.0
113      d_ql(1:ngrid,1:nlayermx) = 0.0
114      d_t(1:ngrid,1:nlayermx)  = 0.0
115      d_q(1:ngrid,1:nlayermx)  = 0.0
[135]116
117!     Calculate v_cptt
118      DO k = 1, nlayermx
[787]119         DO i = 1, ngrid
[135]120            v_cptt(i,k) = RCPD * local_t(i,k)
[728]121            v_t = MAX(local_t(i,k),15.)
[135]122            v_p = pplay(i,k)
123
[728]124            call Psat_water(v_t,v_p,zpsat,zqs(i,k))
125            call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
[135]126         ENDDO
127      ENDDO
128
129!     Calculate Gamma * Cp * dz: (gamma is the critical gradient)
130      DO k = 2, nlayermx
[787]131         DO i = 1, ngrid
[135]132            zdp = pplev(i,k)-pplev(i,k+1)
133            zdpm = pplev(i,k-1)-pplev(i,k)
[728]134            gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
135                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
136                * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
137                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                   
[135]138         ENDDO
139      ENDDO
140
141!------------------------------------ modification of unstable profile
[787]142      DO 9999 i = 1, ngrid
143
[135]144      itest(i) = .FALSE.
145
146!        print*,'we in the loop'
147!        stop   
148
149      k1 = 0
150      k2 = 1
151
152  810 CONTINUE ! look for k1, the base of the column
153      k2 = k2 + 1
154      IF (k2 .GT. nlayermx) GOTO 9999
155      zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
[728]156      zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2))   &
157         +(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
[135]158
159      IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
160      k1 = k2 - 1
161      itest(i) = .TRUE.
162
163  820 CONTINUE !! look for k2, the top of the column
164      IF (k2 .EQ. nlayermx) GOTO 821
165      k2p = k2 + 1
[728]166      zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
[135]167      zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
168
169      IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821
170      k2 = k2p
171      GOTO 820
172  821 CONTINUE
173
174!------------------------------------------------------ local adjustment
175  830 CONTINUE ! actual adjustment
[728]176    Do nn=1,niter
[135]177      v_cptj(k1) = 0.0
178      zdp = pplev(i,k1)-pplev(i,k1+1)
[728]179      v_cptjk1 = ( (1.0+zdqs(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) + RLVTT*(local_q(i,k1)-zqs(i,k1)) ) * zdp
180      v_ssig = zdp * (1.0+zdqs(i,k1))
[135]181
182      k1p = k1 + 1
183      DO k = k1p, k2
184         zdp = pplev(i,k)-pplev(i,k+1)
185         v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
[728]186         v_cptjk1 = v_cptjk1 + zdp * ( (1.0+zdqs(i, k))*(v_cptt(i,k)+v_cptj(k)) + RLVTT*(local_q(i,k)-zqs(i,k)) )       
187         v_ssig = v_ssig + zdp *(1.0+zdqs(i,k))
[135]188      ENDDO
189
190
191      ! this right here is where the adjustment is done???
192      DO k = k1, k2
193         cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
194         cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
[728]195         v_cptt(i,k)=cp_new_t(k)
196         local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT
[135]197         local_t(i,k) = cp_new_t(k) / RCPD
[253]198
[728]199         v_t = MAX(local_t(i,k),15.)
200         v_p = pplay(i,k)
201         
202         call Psat_water(v_t,v_p,zpsat,zqs(i,k))
203         call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
204
205
206
207!          print*,'i,k,zqs,cpdT=',i,k,zqs(i,k),cp_delta_t(k)
[135]208      ENDDO
[728]209    Enddo ! nn=1,niter
[135]210
[253]211
[135]212!--------------------------------------------------- sounding downwards
213!              -- we refine the prognostic variables in
214!              -- the layer about to be adjusted
215
[728]216!      DO k = k1, k2
217!         v_cptt(i,k) = RCPD * local_t(i,k)
218!         v_t = local_t(i,k)
219!         v_p = pplay(i,k)
220!       
221!         call Psat_water(v_t,v_p,zpsat,zqs(i,k))
222!         call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
223!      ENDDO
[135]224
225      DO k = 2, nlayermx
226         zdpm = pplev(i,k-1) - pplev(i,k)
227         zdp = pplev(i,k) - pplev(i,k+1)
[728]228         gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)             &
229                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
230                * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
231                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                   
[135]232      ENDDO
233
234!     Test to see if we've reached the bottom
235
236      IF (k1 .EQ. 1) GOTO 841 ! yes we have!
237      zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
[728]238      zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1))   &
239        + (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
[135]240      IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
241
242  840 CONTINUE
243      k1 = k1 - 1
244      IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
[728]245      zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1))               &
[135]246                  *(pplev(i,k1-1)-pplev(i,k1))
247      zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
248      IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
249         GOTO 840
250      ELSE
251         GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
252      ENDIF
253  841 CONTINUE
254
255      GOTO 810 ! look for other layers higher up
256
257 9999 CONTINUE ! loop over all the points
258
259!      print*,'k1=',k1
260!      print*,'k2=',k2
261
262!      print*,'local_t=',local_t
263!      print*,'v_cptt=',v_cptt
264!      print*,'gamcpdz=',gamcpdz
265
266!-----------------------------------------------------------------------
267! Determine the cloud fraction (hypothese: la nebulosite a lieu
268! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
269
270      DO k = 1, nlayermx
[787]271      DO i = 1, ngrid
[135]272         IF (itest(i)) THEN
[728]273         delta_q(i,k) = local_q(i,k) - zq(i,k)
[135]274         IF (delta_q(i,k).LT.0.) rneb(i,k)  = 1.0
275         ENDIF
276      ENDDO
277      ENDDO
278
279! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
280! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
281! diminue et d'une maniere proportionnelle a cet diminution):
282
[787]283      DO i = 1, ngrid
[135]284         IF (itest(i)) THEN
285         zq1(i) = 0.0
286         zq2(i) = 0.0
287         ENDIF
288      ENDDO
289      DO k = 1, nlayermx
[787]290      DO i = 1, ngrid
[135]291         IF (itest(i)) THEN
292         zdp = pplev(i,k)-pplev(i,k+1)
293         zq1(i) = zq1(i) - delta_q(i,k) * zdp
294         zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
295         ENDIF
296      ENDDO
297      ENDDO
298      DO k = 1, nlayermx
[787]299      DO i = 1, ngrid
[135]300         IF (itest(i)) THEN
[728]301           IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
[135]302         ENDIF
303      ENDDO
304      ENDDO
305
[253]306!      print*,'local_q BEFORE=',local_q
307
[135]308      DO k = 1, nlayermx
[787]309      DO i = 1, ngrid
[135]310          local_q(i, k) = MAX(local_q(i, k), seuil_vap)
311      ENDDO
312      ENDDO
313
314      DO k = 1, nlayermx
[787]315      DO i = 1, ngrid
[728]316         d_t(i,k) = local_t(i,k) - zt(i,k)
317         d_q(i,k) = local_q(i,k) - zq(i,k)
[135]318      ENDDO
319      ENDDO
320
321!     now subroutine -----> GCM variables
322      DO k = 1, nlayermx
[787]323         DO i = 1, ngrid
[135]324           
[728]325            pdtmana(i,k)       = d_t(i,k)/ptimestep
326            pdqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
327            pdqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
[135]328         
329         ENDDO
330      ENDDO
331
[253]332
[135]333      RETURN
[253]334   END
Note: See TracBrowser for help on using the repository browser.