GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: dyn3d/friction.F Lines: 0 35 0.0 %
Date: 2023-06-30 12:56:34 Branches: 0 40 0.0 %

Line Branch Exec Source
1
!
2
! $Id: friction.F 2597 2016-07-22 06:44:47Z emillour $
3
!
4
c=======================================================================
5
      SUBROUTINE friction(ucov,vcov,pdt)
6
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
      USE comconst_mod, ONLY: pi
15
      IMPLICIT 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
      include "dimensions.h"
30
      include "paramet.h"
31
      include "comgeom2.h"
32
      include "iniprint.h"
33
      include "academic.h"
34
35
! arguments:
36
      REAL,INTENT(out) :: ucov( iip1,jjp1,llm )
37
      REAL,INTENT(out) :: vcov( iip1,jjm,llm )
38
      REAL,INTENT(in) :: pdt ! time step
39
40
! local variables:
41
42
      REAL modv(iip1,jjp1),zco,zsi
43
      REAL vpn,vps,upoln,upols,vpols,vpoln
44
      REAL u2(iip1,jjp1),v2(iip1,jjm)
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"
50
      CHARACTER(len=80) :: abort_message
51
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
63
      if (friction_type.eq.0) then
64
c   calcul des composantes au carre du vent naturel
65
      do j=1,jjp1
66
         do i=1,iip1
67
            u2(i,j)=ucov(i,j,1)*ucov(i,j,1)*unscu2(i,j)
68
         enddo
69
      enddo
70
      do j=1,jjm
71
         do i=1,iip1
72
            v2(i,j)=vcov(i,j,1)*vcov(i,j,1)*unscv2(i,j)
73
         enddo
74
      enddo
75
76
c   calcul du module de V en dehors des poles
77
      do j=2,jjm
78
         do i=2,iip1
79
            modv(i,j)=sqrt(0.5*(u2(i-1,j)+u2(i,j)+v2(i,j-1)+v2(i,j)))
80
         enddo
81
         modv(1,j)=modv(iip1,j)
82
      enddo
83
84
c   les deux composantes du vent au pole sont obtenues comme
85
c   premiers modes de fourier de v pres du pole
86
      upoln=0.
87
      vpoln=0.
88
      upols=0.
89
      vpols=0.
90
      do i=2,iip1
91
         zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
92
         zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
93
         vpn=vcov(i,1,1)/cv(i,1)
94
         vps=vcov(i,jjm,1)/cv(i,jjm)
95
         upoln=upoln+zco*vpn
96
         vpoln=vpoln+zsi*vpn
97
         upols=upols+zco*vps
98
         vpols=vpols+zsi*vps
99
      enddo
100
      vpn=sqrt(upoln*upoln+vpoln*vpoln)/pi
101
      vps=sqrt(upols*upols+vpols*vpols)/pi
102
      do i=1,iip1
103
c        modv(i,1)=vpn
104
c        modv(i,jjp1)=vps
105
         modv(i,1)=modv(i,2)
106
         modv(i,jjp1)=modv(i,jjm)
107
      enddo
108
109
c   calcul du frottement au sol.
110
      do j=2,jjm
111
         do i=1,iim
112
            ucov(i,j,1)=ucov(i,j,1)
113
     s      -cfric*pdt*0.5*(modv(i+1,j)+modv(i,j))*ucov(i,j,1)
114
         enddo
115
         ucov(iip1,j,1)=ucov(1,j,1)
116
      enddo
117
      do j=1,jjm
118
         do i=1,iip1
119
            vcov(i,j,1)=vcov(i,j,1)
120
     s      -cfric*pdt*0.5*(modv(i,j+1)+modv(i,j))*vcov(i,j,1)
121
         enddo
122
         vcov(iip1,j,1)=vcov(1,j,1)
123
      enddo
124
      endif ! of if (friction_type.eq.0)
125
126
      if (friction_type.eq.1) then
127
        do l=1,llm
128
          ucov(:,:,l)=ucov(:,:,l)*(1.-pdt*kfrict(l))
129
          vcov(:,:,l)=vcov(:,:,l)*(1.-pdt*kfrict(l))
130
        enddo
131
      endif
132
133
      RETURN
134
      END
135