source: dynamico_lmdz/simple_physics/phyparam/param/convadj.F @ 4183

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

simple_physics : comcstfi.h => MODULE phys_const.F90
File size: 6.2 KB
Line 
1      SUBROUTINE convadj(ngrid,nlay,ptimestep,
2     S                   pplay,pplev,ppopsk,
3     $                   pu,pv,ph,
4     $                   pdufi,pdvfi,pdhfi,
5     $                   pduadj,pdvadj,pdhadj)
6      USE phys_const
7      IMPLICIT NONE
8
9c=======================================================================
10c
11c   ajustement convectif sec
12c   on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement
13c
14c=======================================================================
15
16c-----------------------------------------------------------------------
17c   declarations:
18c   -------------
19
20#include "dimensions.h"
21
22c   arguments:
23c   ----------
24
25      INTEGER ngrid,nlay
26      REAL ptimestep
27      REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
28      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
29      REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
30      REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)
31
32c   local:
33c   ------
34
35      INTEGER ig,i,l,l1,l2,jj
36      INTEGER jcnt, jadrs(ngrid)
37
38      REAL*8 sig(nlay+1),sdsig(nlay),dsig(nlay)
39      REAL*8 zu(ngrid,nlay),zv(ngrid,nlay)
40      REAL*8 zh(ngrid,nlay)
41      REAL*8 zu2(ngrid,nlay),zv2(ngrid,nlay)
42      REAL*8 zh2(ngrid,nlay)
43      REAL*8 zhm,zsm,zum,zvm,zalpha
44
45      LOGICAL vtest(ngrid),down
46
47c
48c-----------------------------------------------------------------------
49c   initialisation:
50c   ---------------
51c
52c
53c-----------------------------------------------------------------------
54c   detection des profils a modifier:
55c   ---------------------------------
56c   si le profil est a modifier
57c   (i.e. ph(niv_sup) < ph(niv_inf) )
58c   alors le tableau "vtest" est mis a .TRUE. ;
59c   sinon, il reste a sa valeur initiale (.FALSE.)
60c   cette operation est vectorisable
61c   On en profite pour copier la valeur initiale de "ph"
62c   dans le champ de travail "zh"
63
64
65      DO 1010 l=1,nlay
66         DO 1015 ig=1,ngrid
67            zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
68            zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
69            zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
701015     CONTINUE
711010  CONTINUE
72
73      zu2(:,:)=zu(:,:)
74      zv2(:,:)=zv(:,:)
75      zh2(:,:)=zh(:,:)
76
77      DO 1020 ig=1,ngrid
78        vtest(ig)=.FALSE.
79 1020 CONTINUE
80c
81      DO 1040 l=2,nlay
82        DO 1060 ig=1,ngrid
83CRAY      vtest(ig)=CVMGM(.TRUE. , vtest(ig),
84CRAY .                      zh2(ig,l)-zh2(ig,l-1))
85          IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE.
86 1060   CONTINUE
87 1040 CONTINUE
88c
89CRAY  CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt)
90      jcnt=0
91      DO 1070 ig=1,ngrid
92         IF(vtest(ig)) THEN
93            jcnt=jcnt+1
94            jadrs(jcnt)=ig
95         ENDIF
961070  CONTINUE
97
98
99c-----------------------------------------------------------------------
100c  Ajustement des "jcnt" profils instables indices par "jadrs":
101c  ------------------------------------------------------------
102c
103      DO 1080 jj = 1, jcnt
104c
105          i = jadrs(jj)
106c
107c   Calcul des niveaux sigma sur cette colonne
108          DO l=1,nlay+1
109            sig(l)=pplev(i,l)/pplev(i,1)
110         ENDDO
111         DO l=1,nlay
112            dsig(l)=sig(l)-sig(l+1)
113            sdsig(l)=ppopsk(i,l)*dsig(l)
114         ENDDO
115          l2 = 1
116c
117c      -- boucle de sondage vers le haut
118c
119cins$     Loop
120 8000     CONTINUE
121c
122            l2 = l2 + 1
123c
124cins$       Exit
125            IF (l2 .GT. nlay) Goto 8001
126c
127            IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN
128c
129c          -- l2 est le niveau le plus haut de la colonne instable
130c
131              l1 = l2 - 1
132              l  = l1
133              zsm = sdsig(l2)
134              zhm = zh2(i, l2)
135c
136c          -- boucle de sondage vers le bas
137c
138cins$         Loop
139 8020         CONTINUE
140c
141                zsm = zsm + sdsig(l)
142                zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
143c
144c            -- doit on etendre la colonne vers le bas ?
145c
146c_EC (M1875) 20/6/87 : AND -> AND THEN
147c
148                down = .FALSE.
149                IF (l1 .NE. 1) THEN    !--  and then
150                  IF (zhm .LT. zh2(i, l1-1)) THEN
151                    down = .TRUE.
152                  END IF
153                END IF
154c
155                IF (down) THEN
156c
157                  l1 = l1 - 1
158                  l  = l1
159c
160                ELSE
161c
162c              -- peut on etendre la colonne vers le haut ?
163c
164cins$             Exit
165                  IF (l2 .EQ. nlay) Goto 8021
166c
167cins$             Exit
168                  IF (zh2(i, l2+1) .GE. zhm) Goto 8021
169c
170                  l2 = l2 + 1
171                  l  = l2
172c
173                END IF
174c
175cins$         End Loop
176              GO TO 8020
177 8021         CONTINUE
178c
179c          -- nouveau profil : constant (valeur moyenne)
180c
181              zalpha=0.
182              zum=0.
183              zvm=0.
184              DO 1100 l = l1, l2
185                zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
186                zh2(i, l) = zhm
187                zum=zum+dsig(l)*zu(i,l)
188                zvm=zvm+dsig(l)*zv(i,l)
189 1100         CONTINUE
190              zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
191              zum=zum/(sig(l1)-sig(l2+1))
192              zvm=zvm/(sig(l1)-sig(l2+1))
193              IF(zalpha.GT.1.) THEN
194                 PRINT*,'WARNING dans convadj zalpha=',zalpha
195                 if(ig.eq.1) then
196                    print*,'Au pole nord'
197                 elseif (ig.eq.ngrid) then
198                    print*,'Au pole sud'
199                 else
200                    print*,'Point i=',
201     .              ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1
202                 endif
203                 STOP
204                 zalpha=1.
205              ELSE
206c                IF(zalpha.LT.0.) STOP'zalpha=0'
207                 IF(zalpha.LT.1.e-5) zalpha=1.e-5
208              ENDIF
209              DO l=l1,l2
210                 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
211                 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
212              ENDDO
213 
214              l2 = l2 + 1
215c
216            END IF
217c
218cins$     End Loop
219          GO TO 8000
220 8001     CONTINUE
221c
222 1080 CONTINUE
223c
224      DO 4000 l=1,nlay
225        DO 4020 ig=1,ngrid
226          pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
227          pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
228          pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
229c         pdhadj(ig,l)=0.
230c         pduadj(ig,l)=0.
231c         pdvadj(ig,l)=0.
232 4020   CONTINUE
233 4000 CONTINUE
234c
235      RETURN
236      END
Note: See TracBrowser for help on using the repository browser.