Changeset 4204 for dynamico_lmdz/simple_physics/phyparam
- Timestamp:
- Dec 20, 2019, 4:14:22 PM (5 years ago)
- Location:
- dynamico_lmdz/simple_physics/phyparam/param
- Files:
-
- 1 edited
- 1 moved
-
convadj.F90 (moved) (moved from dynamico_lmdz/simple_physics/phyparam/param/convadj.F) (1 diff)
-
surface.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
dynamico_lmdz/simple_physics/phyparam/param/convadj.F90
r4196 r4204 1 SUBROUTINE convadj(ngrid,nlay,ptimestep, 2 S pplay,pplev,ppopsk, 3 $ pu,pv,ph, 4 $ pdufi,pdvfi,pdhfi, 5 $ pduadj,pdvadj,pdhadj) 6 USE phys_const 7 IMPLICIT NONE 1 SUBROUTINE convadj(ngrid,nlay,ptimestep, & 2 & pplay,pplev,ppopsk, & 3 & pu,pv,ph, & 4 & pdufi,pdvfi,pdhfi, & 5 & pduadj,pdvadj,pdhadj) 6 USE phys_const 7 IMPLICIT NONE 8 9 !======================================================================= 10 ! 11 ! ajustement convectif sec 12 ! on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement 13 ! 14 !======================================================================= 15 16 !----------------------------------------------------------------------- 17 ! declarations: 18 ! ------------- 19 20 #include "dimensions.h" 21 22 ! arguments: 23 ! ---------- 24 25 INTEGER ngrid,nlay 26 REAL ptimestep 27 REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay) 28 REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay) 29 REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay) 30 REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay) 31 32 ! local: 33 ! ------ 34 35 INTEGER ig,i,l,l1,l2,jj 36 INTEGER jcnt, jadrs(ngrid) 37 38 REAL*8 sig(nlay+1),sdsig(nlay),dsig(nlay) 39 REAL*8 zu(ngrid,nlay),zv(ngrid,nlay) 40 REAL*8 zh(ngrid,nlay) 41 REAL*8 zu2(ngrid,nlay),zv2(ngrid,nlay) 42 REAL*8 zh2(ngrid,nlay) 43 REAL*8 zhm,zsm,zum,zvm,zalpha 44 45 LOGICAL vtest(ngrid),down 46 47 ! 48 !----------------------------------------------------------------------- 49 ! initialisation: 50 ! --------------- 51 ! 52 ! 53 !----------------------------------------------------------------------- 54 ! detection des profils a modifier: 55 ! --------------------------------- 56 ! si le profil est a modifier 57 ! (i.e. ph(niv_sup) < ph(niv_inf) ) 58 ! alors le tableau "vtest" est mis a .TRUE. ; 59 ! sinon, il reste a sa valeur initiale (.FALSE.) 60 ! cette operation est vectorisable 61 ! On en profite pour copier la valeur initiale de "ph" 62 ! dans le champ de travail "zh" 63 64 DO l=1,nlay 65 DO ig=1,ngrid 66 zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep 67 zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep 68 zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep 69 END DO 70 END DO 71 72 zu2(:,:)=zu(:,:) 73 zv2(:,:)=zv(:,:) 74 zh2(:,:)=zh(:,:) 75 76 DO ig=1,ngrid 77 vtest(ig)=.FALSE. 78 END DO 79 ! 80 DO l=2,nlay 81 DO ig=1,ngrid 82 !CRAY vtest(ig)=CVMGM(.TRUE. , vtest(ig), 83 !CRAY . zh2(ig,l)-zh2(ig,l-1)) 84 IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE. 85 END DO 86 END DO 87 ! 88 !CRAY CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt) 89 jcnt=0 90 DO ig=1,ngrid 91 IF(vtest(ig)) THEN 92 jcnt=jcnt+1 93 jadrs(jcnt)=ig 94 ENDIF 95 END DO 96 97 !----------------------------------------------------------------------- 98 ! Ajustement des "jcnt" profils instables indices par "jadrs": 99 ! ------------------------------------------------------------ 100 ! 101 DO jj = 1, jcnt 102 ! 103 i = jadrs(jj) 104 ! 105 ! Calcul des niveaux sigma sur cette colonne 106 DO l=1,nlay+1 107 sig(l)=pplev(i,l)/pplev(i,1) 108 ENDDO 109 DO l=1,nlay 110 dsig(l)=sig(l)-sig(l+1) 111 sdsig(l)=ppopsk(i,l)*dsig(l) 112 ENDDO 113 l2 = 1 114 ! 115 ! -- boucle de sondage vers le haut 116 ! 117 DO WHILE(.TRUE.) 118 ! 8000 CONTINUE 119 ! 120 l2 = l2 + 1 121 ! 122 IF (l2 .GT. nlay) EXIT ! Goto 8001 123 ! 124 IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN 125 ! 126 ! -- l2 est le niveau le plus haut de la colonne instable 127 ! 128 l1 = l2 - 1 129 l = l1 130 zsm = sdsig(l2) 131 zhm = zh2(i, l2) 132 ! 133 ! -- boucle de sondage vers le bas 134 ! 135 DO WHILE(.TRUE.) 136 ! 8020 CONTINUE 137 zsm = zsm + sdsig(l) 138 zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm 139 ! 140 ! -- doit on etendre la colonne vers le bas ? 141 ! 142 !_EC (M1875) 20/6/87 : AND -> AND THEN 143 ! 144 down = .FALSE. 145 IF (l1 .NE. 1) THEN !-- and then 146 IF (zhm .LT. zh2(i, l1-1)) THEN 147 down = .TRUE. 148 END IF 149 END IF 150 ! 151 IF (down) THEN 152 l1 = l1 - 1 153 l = l1 154 ELSE 155 ! -- peut on etendre la colonne vers le haut ? 156 IF (l2 .EQ. nlay) EXIT !Goto 8021 157 IF (zh2(i, l2+1) .GE. zhm) EXIT !Goto 8021 158 l2 = l2 + 1 159 l = l2 160 END IF 161 162 ! GO TO 8020 163 END DO 164 ! 8021 CONTINUE 165 ! 166 ! -- nouveau profil : constant (valeur moyenne) 167 ! 168 zalpha=0. 169 zum=0. 170 zvm=0. 171 DO l = l1, l2 172 zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l) 173 zh2(i, l) = zhm 174 zum=zum+dsig(l)*zu(i,l) 175 zvm=zvm+dsig(l)*zv(i,l) 176 END DO 177 zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1))) 178 zum=zum/(sig(l1)-sig(l2+1)) 179 zvm=zvm/(sig(l1)-sig(l2+1)) 180 IF(zalpha.GT.1.) THEN 181 PRINT*,'WARNING dans convadj zalpha=',zalpha 182 if(ig.eq.1) then 183 print*,'Au pole nord' 184 elseif (ig.eq.ngrid) then 185 print*,'Au pole sud' 186 else 187 print*,'Point i=', & 188 ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1 189 endif 190 STOP 191 zalpha=1. 192 ELSE 193 ! IF(zalpha.LT.0.) STOP'zalpha=0' 194 IF(zalpha.LT.1.e-5) zalpha=1.e-5 195 ENDIF 196 DO l=l1,l2 197 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l)) 198 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) 199 ENDDO 200 201 l2 = l2 + 1 202 ! 203 END IF 204 ! 205 ! GO TO 8000 206 END DO 207 ! 8001 CONTINUE 208 ! 209 ! 210 DO l=1,nlay 211 DO ig=1,ngrid 212 pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep 213 pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep 214 pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep 215 ! pdhadj(ig,l)=0. 216 ! pduadj(ig,l)=0. 217 ! pdvadj(ig,l)=0. 218 END DO 219 END DO 220 ! 221 END DO 8 222 9 c======================================================================= 10 c 11 c ajustement convectif sec 12 c on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement 13 c 14 c======================================================================= 15 16 c----------------------------------------------------------------------- 17 c declarations: 18 c ------------- 19 20 #include "dimensions.h" 21 22 c arguments: 23 c ---------- 24 25 INTEGER ngrid,nlay 26 REAL ptimestep 27 REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay) 28 REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay) 29 REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay) 30 REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay) 31 32 c local: 33 c ------ 34 35 INTEGER ig,i,l,l1,l2,jj 36 INTEGER jcnt, jadrs(ngrid) 37 38 REAL*8 sig(nlay+1),sdsig(nlay),dsig(nlay) 39 REAL*8 zu(ngrid,nlay),zv(ngrid,nlay) 40 REAL*8 zh(ngrid,nlay) 41 REAL*8 zu2(ngrid,nlay),zv2(ngrid,nlay) 42 REAL*8 zh2(ngrid,nlay) 43 REAL*8 zhm,zsm,zum,zvm,zalpha 44 45 LOGICAL vtest(ngrid),down 46 47 c 48 c----------------------------------------------------------------------- 49 c initialisation: 50 c --------------- 51 c 52 c 53 c----------------------------------------------------------------------- 54 c detection des profils a modifier: 55 c --------------------------------- 56 c si le profil est a modifier 57 c (i.e. ph(niv_sup) < ph(niv_inf) ) 58 c alors le tableau "vtest" est mis a .TRUE. ; 59 c sinon, il reste a sa valeur initiale (.FALSE.) 60 c cette operation est vectorisable 61 c On en profite pour copier la valeur initiale de "ph" 62 c dans le champ de travail "zh" 63 64 65 DO 1010 l=1,nlay 66 DO 1015 ig=1,ngrid 67 zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep 68 zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep 69 zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep 70 1015 CONTINUE 71 1010 CONTINUE 72 73 zu2(:,:)=zu(:,:) 74 zv2(:,:)=zv(:,:) 75 zh2(:,:)=zh(:,:) 76 77 DO 1020 ig=1,ngrid 78 vtest(ig)=.FALSE. 79 1020 CONTINUE 80 c 81 DO 1040 l=2,nlay 82 DO 1060 ig=1,ngrid 83 CRAY vtest(ig)=CVMGM(.TRUE. , vtest(ig), 84 CRAY . zh2(ig,l)-zh2(ig,l-1)) 85 IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE. 86 1060 CONTINUE 87 1040 CONTINUE 88 c 89 CRAY CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt) 90 jcnt=0 91 DO 1070 ig=1,ngrid 92 IF(vtest(ig)) THEN 93 jcnt=jcnt+1 94 jadrs(jcnt)=ig 95 ENDIF 96 1070 CONTINUE 97 98 99 c----------------------------------------------------------------------- 100 c Ajustement des "jcnt" profils instables indices par "jadrs": 101 c ------------------------------------------------------------ 102 c 103 DO 1080 jj = 1, jcnt 104 c 105 i = jadrs(jj) 106 c 107 c Calcul des niveaux sigma sur cette colonne 108 DO l=1,nlay+1 109 sig(l)=pplev(i,l)/pplev(i,1) 110 ENDDO 111 DO l=1,nlay 112 dsig(l)=sig(l)-sig(l+1) 113 sdsig(l)=ppopsk(i,l)*dsig(l) 114 ENDDO 115 l2 = 1 116 c 117 c -- boucle de sondage vers le haut 118 c 119 cins$ Loop 120 8000 CONTINUE 121 c 122 l2 = l2 + 1 123 c 124 cins$ Exit 125 IF (l2 .GT. nlay) Goto 8001 126 c 127 IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN 128 c 129 c -- l2 est le niveau le plus haut de la colonne instable 130 c 131 l1 = l2 - 1 132 l = l1 133 zsm = sdsig(l2) 134 zhm = zh2(i, l2) 135 c 136 c -- boucle de sondage vers le bas 137 c 138 cins$ Loop 139 8020 CONTINUE 140 c 141 zsm = zsm + sdsig(l) 142 zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm 143 c 144 c -- doit on etendre la colonne vers le bas ? 145 c 146 c_EC (M1875) 20/6/87 : AND -> AND THEN 147 c 148 down = .FALSE. 149 IF (l1 .NE. 1) THEN !-- and then 150 IF (zhm .LT. zh2(i, l1-1)) THEN 151 down = .TRUE. 152 END IF 153 END IF 154 c 155 IF (down) THEN 156 c 157 l1 = l1 - 1 158 l = l1 159 c 160 ELSE 161 c 162 c -- peut on etendre la colonne vers le haut ? 163 c 164 cins$ Exit 165 IF (l2 .EQ. nlay) Goto 8021 166 c 167 cins$ Exit 168 IF (zh2(i, l2+1) .GE. zhm) Goto 8021 169 c 170 l2 = l2 + 1 171 l = l2 172 c 173 END IF 174 c 175 cins$ End Loop 176 GO TO 8020 177 8021 CONTINUE 178 c 179 c -- nouveau profil : constant (valeur moyenne) 180 c 181 zalpha=0. 182 zum=0. 183 zvm=0. 184 DO 1100 l = l1, l2 185 zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l) 186 zh2(i, l) = zhm 187 zum=zum+dsig(l)*zu(i,l) 188 zvm=zvm+dsig(l)*zv(i,l) 189 1100 CONTINUE 190 zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1))) 191 zum=zum/(sig(l1)-sig(l2+1)) 192 zvm=zvm/(sig(l1)-sig(l2+1)) 193 IF(zalpha.GT.1.) THEN 194 PRINT*,'WARNING dans convadj zalpha=',zalpha 195 if(ig.eq.1) then 196 print*,'Au pole nord' 197 elseif (ig.eq.ngrid) then 198 print*,'Au pole sud' 199 else 200 print*,'Point i=', 201 . ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1 202 endif 203 STOP 204 zalpha=1. 205 ELSE 206 c IF(zalpha.LT.0.) STOP'zalpha=0' 207 IF(zalpha.LT.1.e-5) zalpha=1.e-5 208 ENDIF 209 DO l=l1,l2 210 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l)) 211 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) 212 ENDDO 213 214 l2 = l2 + 1 215 c 216 END IF 217 c 218 cins$ End Loop 219 GO TO 8000 220 8001 CONTINUE 221 c 222 1080 CONTINUE 223 c 224 DO 4000 l=1,nlay 225 DO 4020 ig=1,ngrid 226 pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep 227 pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep 228 pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep 229 c pdhadj(ig,l)=0. 230 c pduadj(ig,l)=0. 231 c pdvadj(ig,l)=0. 232 4020 CONTINUE 233 4000 CONTINUE 234 c 235 RETURN 236 END 223 END SUBROUTINE convadj -
dynamico_lmdz/simple_physics/phyparam/param/surface.h
r4176 r4204 1 ! REAL rnaturfi(ngridmax),phisfi(ngridmax)2 1 REAL I_mer,I_ter,Cd_mer,Cd_ter 3 2 REAL alb_mer,alb_ter,emi_mer,emi_ter 4 3 5 4 COMMON/surface/ 6 ! rnaturfi,phisfi,7 5 s I_mer,I_ter,Cd_mer,Cd_ter, 8 6 s alb_mer,alb_ter,emi_mer,emi_ter
Note: See TracChangeset
for help on using the changeset viewer.
