source: dynamico_lmdz/simple_physics/phyparam/param/coefdifv.F @ 4176

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

simple_physics : copy code from FH

File size: 2.2 KB
Line 
1      SUBROUTINE coefdifv( np,nlev,
2     $               pu,pv,ph,pphi,
3     $               pcdzv,pcdzh)
4      USE constlim
5      IMPLICIT NONE
6
7c=======================================================================
8c
9c   calcul des coeficients de la diffusion verticale
10c
11c=======================================================================
12c
13c-----------------------------------------------------------------------
14c   Declarations:
15c   -------------
16
17#include "dimensions.h"
18#include "comcstfi.h"
19
20c   Arguments:
21c   ----------
22
23      INTEGER np,nlev
24      REAL pu(np,nlev),pv(np,nlev)
25      REAL ph(np,nlev),pphi(np,nlev)
26      REAL pcdzv(np,nlev),pcdzh(np,nlev)
27
28c   Local:
29c   ------
30
31      INTEGER ilev,ip
32      REAL z1(np),z2(np)
33      REAL z1odz(np),zdzu2(np),zhbar(np)
34
35c-----------------------------------------------------------------------
36c   initialisations:
37c   ----------------
38
39c-----------------------------------------------------------------------
40c   couche de surface:
41c   ------------------
42
43      DO 210 ip=1,np
44         z1(ip)=pu(ip,1)*pu(ip,1)+pv(ip,1)*pv(ip,1)
45         pcdzv(ip,1)=dgcdrag(ip)*(SQRT(z1(ip))+1.)  / ph(ip,1)
46         pcdzh(ip,1)=pcdzv(ip,1)
47210   CONTINUE
48c
49c
50c-----------------------------------------------------------------------
51c   autres couches:
52c   ---------------
53
54      PRINT*,'Coeff k'
55      DO 310 ilev=1,nlev-1
56         DO 320 ip=1,np
57            z1(ip)=pu(ip,ilev+1)-pu(ip,ilev)
58            z2(ip)=pv(ip,ilev+1)-pv(ip,ilev)
59            z1(ip)=z1(ip)*z1(ip)+z2(ip)*z2(ip)
60            z1odz(ip)=g/(pphi(ip,ilev+1)-pphi(ip,ilev))
61            zdzu2(ip)=z1(ip)*z1odz(ip)*z1odz(ip)
62            zhbar(ip)=0.5*(ph(ip,ilev)+ph(ip,ilev+1))
63            z1(ip)=( (ph(ip,ilev+1)-ph(ip,ilev))*z1odz(ip)-
64     $             cpgam ) * ccdzh/zhbar(ip)
65            pcdzv(ip,ilev+1)=cdzconst(ilev+1)*
66     $            SQRT(AMAX1(zdzu2(ip)-z1(ip),cdzmin)) /
67c    $            SQRT(cdzmin) /
68     $            (zhbar(ip)*zhbar(ip))
69            IF(ip.eq.np/2+1) PRINT*,lmixmin*lmixmin*
70     s     SQRT(AMAX1(zdzu2(ip)-z1(ip),cdzmin))
71            pcdzh(ip,ilev+1)=pcdzv(ip,ilev+1)
72320      CONTINUE
73310   CONTINUE
74
75c-----------------------------------------------------------------------
76
77      RETURN
78      END
Note: See TracBrowser for help on using the repository browser.