Changeset 5246 for LMDZ6/trunk/libf/dyn3dmem/friction_loc.F90
- Timestamp:
- Oct 21, 2024, 2:58:45 PM (23 hours ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3dmem/friction_loc.F90
r5245 r5246 2 2 ! $Id: friction_p.F 1299 2010-01-20 14:27:21Z fairhead $ 3 3 ! 4 c=======================================================================5 6 7 4 !======================================================================= 5 SUBROUTINE friction_loc(ucov,vcov,pdt) 6 USE parallel_lmdz 7 USE control_mod 8 8 #ifdef CPP_IOIPSL 9 9 USE IOIPSL 10 10 #else 11 ! if not using IOIPSL, we still need to use (a local version of) getin12 11 ! if not using IOIPSL, we still need to use (a local version of) getin 12 USE ioipsl_getincom 13 13 #endif 14 15 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 lowermost23 !atmospheric layer24 !friction_type=1 : Friction applied on all atmospheric layer (but25 !(default) with stronger magnitude near the surface; see26 !iniacademic.F)27 !=======================================================================28 29 30 31 32 33 34 35 ! arguments:36 37 38 39 40 ! local variables:41 42 REALmodv(iip1,jjb_u:jje_u),zco,zsi43 REALvpn,vps,upoln,upols,vpols,vpoln44 REALu2(iip1,jjb_u:jje_u),v2(iip1,jjb_v:jje_v)45 INTEGERi,j,l46 47 48 49 50 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(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 51 !$OMP THREADPRIVATE(firstcall,friction_type) 52 52 integer :: jjb,jje 53 53 54 54 !$OMP SINGLE 55 56 57 58 59 60 61 62 63 64 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 65 !$OMP END SINGLE COPYPRIVATE(friction_type,firstcall) 66 66 67 67 if (friction_type.eq.0) then ! friction on first layer only 68 68 !$OMP SINGLE 69 ccalcul des composantes au carre du vent naturel70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 ccalcul du module de V en dehors des poles92 93 94 95 96 97 98 99 100 101 102 103 104 cles deux composantes du vent au pole sont obtenues comme105 cpremiers modes de fourier de v pres du pole106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 cmodv(i,1)=vpn121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 cmodv(i,jjp1)=vps140 141 142 143 144 145 ccalcul du frottement au sol.146 147 148 149 150 151 152 153 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 157 158 159 160 161 162 163 164 165 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 169 170 69 ! 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 91 ! 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 104 ! les deux composantes du vent au pole sont obtenues comme 105 ! 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 120 ! 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 139 ! modv(i,jjp1)=vps 140 modv(i,jjp1)=modv(i,jjm) 141 enddo 142 143 endif 144 145 ! 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 -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 -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 171 !$OMP END SINGLE 172 173 174 175 ! for ucov()176 177 178 179 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 180 181 181 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 182 183 ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)*184 &(1.-pdt*kfrict(l))185 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 186 !$OMP END DO NOWAIT 187 188 189 190 191 192 187 188 ! ! for vcoc() 189 jjb=jj_begin 190 jje=jj_end 191 if (pole_sud) jje=jj_end-1 192 193 193 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 194 195 vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)*196 &(1.-pdt*kfrict(l))197 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 198 !$OMP END DO 199 200 201 202 END 203 199 endif ! of if (friction_type.eq.1) 200 201 RETURN 202 END SUBROUTINE friction_loc 203
Note: See TracChangeset
for help on using the changeset viewer.