source: LMDZ5/branches/testing/libf/dyn3dmem/friction_loc.F @ 1707

Last change on this file since 1707 was 1707, checked in by Laurent Fairhead, 11 years ago

Version testing basée sur la r1706

Testing release based on r1706
File size: 5.4 KB
Line 
1!
2! $Id: friction_p.F 1299 2010-01-20 14:27:21Z fairhead $
3!
4c=======================================================================
5      SUBROUTINE friction_loc(ucov,vcov,pdt)
6      USE parallel
7      USE control_mod
8#ifdef CPP_IOIPSL
9      USE IOIPSL
10#else
11! if not using IOIPSL, we still need to use (a local version of) getin
12      USE ioipsl_getincom
13#endif
14      IMPLICIT NONE
15
16!=======================================================================
17!
18!   Friction for the Newtonian case:
19!   --------------------------------
20!    2 possibilities (depending on flag 'friction_type'
21!     friction_type=0 : A friction that is only applied to the lowermost
22!                       atmospheric layer
23!     friction_type=1 : Friction applied on all atmospheric layer (but
24!       (default)       with stronger magnitude near the surface; see
25!                       iniacademic.F)
26!=======================================================================
27
28#include "dimensions.h"
29#include "paramet.h"
30#include "comgeom2.h"
31#include "comconst.h"
32#include "iniprint.h"
33#include "academic.h"
34
35! arguments:
36      REAL,INTENT(inout) :: ucov( iip1,jjb_u:jje_u,llm )
37      REAL,INTENT(inout) :: vcov( iip1,jjb_v:jje_v,llm )
38      REAL,INTENT(in) :: pdt ! time step
39
40! local variables:
41
42      REAL modv(iip1,jjb_u:jje_u),zco,zsi
43      REAL vpn,vps,upoln,upols,vpols,vpoln
44      REAL u2(iip1,jjb_u:jje_u),v2(iip1,jjb_v:jje_v)
45      INTEGER  i,j,l
46      REAL,PARAMETER :: cfric=1.e-5
47      LOGICAL,SAVE :: firstcall=.true.
48      INTEGER,SAVE :: friction_type=1
49      CHARACTER(len=20) :: modname="friction_p"
50      CHARACTER(len=80) :: abort_message
51!$OMP THREADPRIVATE(firstcall,friction_type)
52      integer :: jjb,jje
53
54!$OMP SINGLE
55      IF (firstcall) THEN
56        ! set friction type
57        call getin("friction_type",friction_type)
58        if ((friction_type.lt.0).or.(friction_type.gt.1)) then
59          abort_message="wrong friction type"
60          write(lunout,*)'Friction: wrong friction type',friction_type
61          call abort_gcm(modname,abort_message,42)
62        endif
63        firstcall=.false.
64      ENDIF
65!$OMP END SINGLE COPYPRIVATE(friction_type,firstcall)
66
67      if (friction_type.eq.0) then ! friction on first layer only
68!$OMP SINGLE
69c   calcul des composantes au carre du vent naturel
70      jjb=jj_begin
71      jje=jj_end+1
72      if (pole_sud) jje=jj_end
73     
74      do j=jjb,jje
75         do i=1,iip1
76            u2(i,j)=ucov(i,j,1)*ucov(i,j,1)*unscu2(i,j)
77         enddo
78      enddo
79     
80      jjb=jj_begin-1
81      jje=jj_end+1
82      if (pole_nord) jjb=jj_begin
83      if (pole_sud) jje=jj_end-1
84     
85      do j=jjb,jje
86         do i=1,iip1
87            v2(i,j)=vcov(i,j,1)*vcov(i,j,1)*unscv2(i,j)
88         enddo
89      enddo
90
91c   calcul du module de V en dehors des poles
92      jjb=jj_begin
93      jje=jj_end+1
94      if (pole_nord) jjb=jj_begin+1
95      if (pole_sud) jje=jj_end-1
96     
97      do j=jjb,jje
98         do i=2,iip1
99            modv(i,j)=sqrt(0.5*(u2(i-1,j)+u2(i,j)+v2(i,j-1)+v2(i,j)))
100         enddo
101         modv(1,j)=modv(iip1,j)
102      enddo
103
104c   les deux composantes du vent au pole sont obtenues comme
105c   premiers modes de fourier de v pres du pole
106      if (pole_nord) then
107     
108        upoln=0.
109        vpoln=0.
110     
111        do i=2,iip1
112           zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
113           zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
114           vpn=vcov(i,1,1)/cv(i,1)
115           upoln=upoln+zco*vpn
116           vpoln=vpoln+zsi*vpn
117        enddo
118        vpn=sqrt(upoln*upoln+vpoln*vpoln)/pi
119        do i=1,iip1
120c          modv(i,1)=vpn
121           modv(i,1)=modv(i,2)
122        enddo
123
124      endif
125     
126      if (pole_sud) then
127     
128        upols=0.
129        vpols=0.
130        do i=2,iip1
131           zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
132           zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
133           vps=vcov(i,jjm,1)/cv(i,jjm)
134           upols=upols+zco*vps
135           vpols=vpols+zsi*vps
136        enddo
137        vps=sqrt(upols*upols+vpols*vpols)/pi
138        do i=1,iip1
139c        modv(i,jjp1)=vps
140         modv(i,jjp1)=modv(i,jjm)
141        enddo
142     
143      endif
144     
145c   calcul du frottement au sol.
146
147      jjb=jj_begin
148      jje=jj_end
149      if (pole_nord) jjb=jj_begin+1
150      if (pole_sud) jje=jj_end-1
151
152      do j=jjb,jje
153         do i=1,iim
154            ucov(i,j,1)=ucov(i,j,1)
155     s      -cfric*pdt*0.5*(modv(i+1,j)+modv(i,j))*ucov(i,j,1)
156         enddo
157         ucov(iip1,j,1)=ucov(1,j,1)
158      enddo
159     
160      jjb=jj_begin
161      jje=jj_end
162      if (pole_sud) jje=jj_end-1
163     
164      do j=jjb,jje
165         do i=1,iip1
166            vcov(i,j,1)=vcov(i,j,1)
167     s      -cfric*pdt*0.5*(modv(i,j+1)+modv(i,j))*vcov(i,j,1)
168         enddo
169         vcov(iip1,j,1)=vcov(1,j,1)
170      enddo
171!$OMP END SINGLE
172      endif ! of if (friction_type.eq.0)
173
174      if (friction_type.eq.1) then
175       ! for ucov()
176        jjb=jj_begin
177        jje=jj_end
178        if (pole_nord) jjb=jj_begin+1
179        if (pole_sud) jje=jj_end-1
180
181!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
182        do l=1,llm
183          ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)*
184     &                                  (1.-pdt*kfrict(l))
185        enddo
186!$OMP END DO NOWAIT
187       
188       ! for vcoc()
189        jjb=jj_begin
190        jje=jj_end
191        if (pole_sud) jje=jj_end-1
192       
193!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
194        do l=1,llm
195          vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)*
196     &                                  (1.-pdt*kfrict(l))
197        enddo
198!$OMP END DO
199      endif ! of if (friction_type.eq.1)
200
201      RETURN
202      END
203
Note: See TracBrowser for help on using the repository browser.