source: LMDZ5/trunk/libf/dyn3dpar/friction_p.F @ 2239

Last change on this file since 2239 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.3 KB
Line 
1!
2! $Id: friction_p.F 1907 2013-11-26 13:10:46Z emillour $
3!
4c=======================================================================
5      SUBROUTINE friction_p(ucov,vcov,pdt)
6      USE parallel_lmdz
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,jjp1,llm )
37      REAL,INTENT(inout) :: vcov( iip1,jjm,llm )
38      REAL,INTENT(in) :: pdt ! time step
39
40! local variables:
41      REAL modv(iip1,jjp1),zco,zsi
42      REAL vpn,vps,upoln,upols,vpols,vpoln
43      REAL u2(iip1,jjp1),v2(iip1,jjm)
44      INTEGER  i,j,l
45      REAL,PARAMETER :: cfric=1.e-5
46      LOGICAL,SAVE :: firstcall=.true.
47      INTEGER,SAVE :: friction_type=1
48      CHARACTER(len=20) :: modname="friction_p"
49      CHARACTER(len=80) :: abort_message
50!$OMP THREADPRIVATE(firstcall,friction_type)
51      integer :: jjb,jje
52
53!$OMP SINGLE
54      IF (firstcall) THEN
55        ! set friction type
56        call getin("friction_type",friction_type)
57        if ((friction_type.lt.0).or.(friction_type.gt.1)) then
58          abort_message="wrong friction type"
59          write(lunout,*)'Friction: wrong friction type',friction_type
60          call abort_gcm(modname,abort_message,42)
61        endif
62        firstcall=.false.
63      ENDIF
64!$OMP END SINGLE COPYPRIVATE(friction_type,firstcall)
65
66      if (friction_type.eq.0) then ! friction on first layer only
67!$OMP SINGLE
68c   calcul des composantes au carre du vent naturel
69      jjb=jj_begin
70      jje=jj_end+1
71      if (pole_sud) jje=jj_end
72     
73      do j=jjb,jje
74         do i=1,iip1
75            u2(i,j)=ucov(i,j,1)*ucov(i,j,1)*unscu2(i,j)
76         enddo
77      enddo
78     
79      jjb=jj_begin-1
80      jje=jj_end+1
81      if (pole_nord) jjb=jj_begin
82      if (pole_sud) jje=jj_end-1
83     
84      do j=jjb,jje
85         do i=1,iip1
86            v2(i,j)=vcov(i,j,1)*vcov(i,j,1)*unscv2(i,j)
87         enddo
88      enddo
89
90c   calcul du module de V en dehors des poles
91      jjb=jj_begin
92      jje=jj_end+1
93      if (pole_nord) jjb=jj_begin+1
94      if (pole_sud) jje=jj_end-1
95     
96      do j=jjb,jje
97         do i=2,iip1
98            modv(i,j)=sqrt(0.5*(u2(i-1,j)+u2(i,j)+v2(i,j-1)+v2(i,j)))
99         enddo
100         modv(1,j)=modv(iip1,j)
101      enddo
102
103c   les deux composantes du vent au pole sont obtenues comme
104c   premiers modes de fourier de v pres du pole
105      if (pole_nord) then
106     
107        upoln=0.
108        vpoln=0.
109     
110        do i=2,iip1
111           zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
112           zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
113           vpn=vcov(i,1,1)/cv(i,1)
114           upoln=upoln+zco*vpn
115           vpoln=vpoln+zsi*vpn
116        enddo
117        vpn=sqrt(upoln*upoln+vpoln*vpoln)/pi
118        do i=1,iip1
119c          modv(i,1)=vpn
120           modv(i,1)=modv(i,2)
121        enddo
122
123      endif
124     
125      if (pole_sud) then
126     
127        upols=0.
128        vpols=0.
129        do i=2,iip1
130           zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
131           zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
132           vps=vcov(i,jjm,1)/cv(i,jjm)
133           upols=upols+zco*vps
134           vpols=vpols+zsi*vps
135        enddo
136        vps=sqrt(upols*upols+vpols*vpols)/pi
137        do i=1,iip1
138c        modv(i,jjp1)=vps
139         modv(i,jjp1)=modv(i,jjm)
140        enddo
141     
142      endif
143     
144c   calcul du frottement au sol.
145
146      jjb=jj_begin
147      jje=jj_end
148      if (pole_nord) jjb=jj_begin+1
149      if (pole_sud) jje=jj_end-1
150
151      do j=jjb,jje
152         do i=1,iim
153            ucov(i,j,1)=ucov(i,j,1)
154     s      -cfric*pdt*0.5*(modv(i+1,j)+modv(i,j))*ucov(i,j,1)
155         enddo
156         ucov(iip1,j,1)=ucov(1,j,1)
157      enddo
158     
159      jjb=jj_begin
160      jje=jj_end
161      if (pole_sud) jje=jj_end-1
162     
163      do j=jjb,jje
164         do i=1,iip1
165            vcov(i,j,1)=vcov(i,j,1)
166     s      -cfric*pdt*0.5*(modv(i,j+1)+modv(i,j))*vcov(i,j,1)
167         enddo
168         vcov(iip1,j,1)=vcov(1,j,1)
169      enddo
170!$OMP END SINGLE
171      endif ! of if (friction_type.eq.0)
172
173      if (friction_type.eq.1) then
174       ! for ucov()
175        jjb=jj_begin
176        jje=jj_end
177        if (pole_nord) jjb=jj_begin+1
178        if (pole_sud) jje=jj_end-1
179
180!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
181        do l=1,llm
182          ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)*
183     &                                  (1.-pdt*kfrict(l))
184        enddo
185!$OMP END DO NOWAIT
186       
187       ! for vcoc()
188        jjb=jj_begin
189        jje=jj_end
190        if (pole_sud) jje=jj_end-1
191       
192!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
193        do l=1,llm
194          vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)*
195     &                                  (1.-pdt*kfrict(l))
196        enddo
197!$OMP END DO
198      endif ! of if (friction_type.eq.1)
199
200      RETURN
201      END
202
Note: See TracBrowser for help on using the repository browser.