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

Last change on this file since 5292 was 5292, checked in by abarral, 4 days ago

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