source: LMDZ6/trunk/libf/dyn3dmem/friction_loc.f90 @ 5280

Last change on this file since 5280 was 5272, checked in by abarral, 7 weeks ago

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