Changeset 1489 for LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d
- Timestamp:
- Feb 17, 2011, 7:20:44 PM (14 years ago)
- Location:
- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d
- Files:
-
- 1 edited
- 1 moved
-
comdissipn.h (modified) (1 diff)
-
inidissip.F90 (moved) (moved from LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/inidissip.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/comdissipn.h
r524 r1489 2 2 ! $Header$ 3 3 ! 4 c----------------------------------------------------------------------- 5 c INCLUDE comdissipn.h 4 ! Attention : ce fichier include est compatible format fixe/format libre 5 ! veillez à n'utiliser que des ! pour les commentaires 6 ! et à bien positionner les & des lignes de continuation 7 ! (les placer en colonne 6 et en colonne 73) 8 !----------------------------------------------------------------------- 9 ! INCLUDE comdissipn.h 6 10 7 11 REAL tetaudiv, tetaurot, tetah, cdivu, crot, cdivh 8 c 9 COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm) , 10 1cdivu, crot, cdivh12 ! 13 COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm) , & 14 & cdivu, crot, cdivh 11 15 12 c 13 cLes parametres de ce common proviennent des calculs effectues dans14 cInidissip .15 c 16 c-----------------------------------------------------------------------16 ! 17 ! Les parametres de ce common proviennent des calculs effectues dans 18 ! Inidissip . 19 ! 20 !----------------------------------------------------------------------- -
LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/inidissip.F90
r1488 r1489 2 2 ! $Id$ 3 3 ! 4 SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , 5 *tetagdiv,tetagrot,tetatemp )6 c=======================================================================7 cinitialisation de la dissipation horizontale8 c=======================================================================9 c-----------------------------------------------------------------------10 cdeclarations:11 c-------------12 13 14 15 16 #include "dimensions.h"17 #include "paramet.h"18 #include "comdissipn.h"19 #include "comconst.h"20 #include "comvert.h"21 #include "logic.h"22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 c-----------------------------------------------------------------------39 c 40 ccalcul des valeurs propres des operateurs par methode iterrative:41 c-----------------------------------------------------------------42 43 44 45 46 47 ccalcul de la valeur propre de divgrad:48 c--------------------------------------49 50 51 4 SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , & 5 tetagdiv,tetagrot,tetatemp ) 6 !======================================================================= 7 ! initialisation de la dissipation horizontale 8 !======================================================================= 9 !----------------------------------------------------------------------- 10 ! declarations: 11 ! ------------- 12 13 USE control_mod 14 15 IMPLICIT NONE 16 include "dimensions.h" 17 include "paramet.h" 18 include "comdissipn.h" 19 include "comconst.h" 20 include "comvert.h" 21 include "logic.h" 22 23 LOGICAL lstardis 24 INTEGER nitergdiv,nitergrot,niterh 25 REAL tetagdiv,tetagrot,tetatemp 26 REAL fact,zvert(llm),zz 27 REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm) 28 REAL ullm,vllm,umin,vmin,zhmin,zhmax 29 REAL zllm,z1llm 30 31 INTEGER l,ij,idum,ii 32 REAL tetamin 33 REAL pseudoz 34 35 REAL ran1 36 37 38 !----------------------------------------------------------------------- 39 ! 40 ! calcul des valeurs propres des operateurs par methode iterrative: 41 ! ----------------------------------------------------------------- 42 43 crot = -1. 44 cdivu = -1. 45 cdivh = -1. 46 47 ! calcul de la valeur propre de divgrad: 48 ! -------------------------------------- 49 idum = 0 50 DO l = 1, llm 51 DO ij = 1, ip1jmp1 52 52 deltap(ij,l) = 1. 53 ENDDO 54 ENDDO 55 56 idum = -1 57 zh(1) = RAN1(idum)-.5 58 idum = 0 59 DO ij = 2, ip1jmp1 60 zh(ij) = RAN1(idum) -.5 61 ENDDO 62 63 CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) 64 65 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 66 67 IF ( zhmin .GE. zhmax ) THEN 68 PRINT*,' Inidissip zh min max ',zhmin,zhmax 69 STOP'probleme generateur alleatoire dans inidissip' 70 ENDIF 71 72 zllm = ABS( zhmax ) 73 DO l = 1,50 74 IF(lstardis) THEN 75 CALL divgrad2(1,zh,deltap,niterh,zh) 76 ELSE 77 CALL divgrad (1,zh,niterh,zh) 78 ENDIF 79 80 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 81 82 zllm = ABS( zhmax ) 83 z1llm = 1./zllm 84 DO ij = 1,ip1jmp1 85 zh(ij) = zh(ij)* z1llm 86 ENDDO 87 ENDDO 88 89 IF(lstardis) THEN 90 cdivh = 1./ zllm 91 ELSE 92 cdivh = zllm ** ( -1./niterh ) 93 ENDIF 94 95 c calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) 96 c ----------------------------------------------------------------- 97 print*,'calcul des valeurs propres' 98 99 DO 20 ii = 1, 2 100 c 101 DO ij = 1, ip1jmp1 102 zu(ij) = RAN1(idum) -.5 103 ENDDO 104 CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) 105 DO ij = 1, ip1jm 106 zv(ij) = RAN1(idum) -.5 107 ENDDO 108 CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) 109 110 CALL minmax(iip1*jjp1,zu,umin,ullm ) 111 CALL minmax(iip1*jjm, zv,vmin,vllm ) 112 113 ullm = ABS ( ullm ) 114 vllm = ABS ( vllm ) 115 116 DO 5 l = 1, 50 117 IF(ii.EQ.1) THEN 118 ccccc CALL covcont( 1,zu,zv,zu,zv ) 119 IF(lstardis) THEN 120 CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv ) 121 ELSE 122 CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv ) 123 ENDIF 124 ELSE 125 IF(lstardis) THEN 126 CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv ) 127 ELSE 128 CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv ) 129 ENDIF 130 ENDIF 131 132 CALL minmax(iip1*jjp1,zu,umin,ullm ) 133 CALL minmax(iip1*jjm, zv,vmin,vllm ) 134 135 ullm = ABS ( ullm ) 136 vllm = ABS ( vllm ) 137 138 zllm = MAX( ullm,vllm ) 139 z1llm = 1./ zllm 140 DO ij = 1, ip1jmp1 141 zu(ij) = zu(ij)* z1llm 142 ENDDO 143 DO ij = 1, ip1jm 144 zv(ij) = zv(ij)* z1llm 145 ENDDO 146 5 CONTINUE 147 148 IF ( ii.EQ.1 ) THEN 149 IF(lstardis) THEN 150 cdivu = 1./zllm 151 ELSE 152 cdivu = zllm **( -1./nitergdiv ) 153 ENDIF 154 ELSE 155 IF(lstardis) THEN 156 crot = 1./ zllm 157 ELSE 158 crot = zllm **( -1./nitergrot ) 159 ENDIF 160 ENDIF 161 162 20 CONTINUE 163 164 c petit test pour les operateurs non star: 165 c ---------------------------------------- 166 167 c IF(.NOT.lstardis) THEN 168 fact = rad*24./REAL(jjm) 169 fact = fact*fact 170 PRINT*,'coef u ', fact/cdivu, 1./cdivu 171 PRINT*,'coef r ', fact/crot , 1./crot 172 PRINT*,'coef h ', fact/cdivh, 1./cdivh 173 c ENDIF 174 175 c----------------------------------------------------------------------- 176 c variation verticale du coefficient de dissipation: 177 c -------------------------------------------------- 178 179 if (ok_strato .and. llm==39) then 180 do l=1,llm 181 pseudoz=8.*log(preff/presnivs(l)) 182 zvert(l)=1+ 183 s (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. 184 s *(dissip_factz-1.) 185 enddo 186 else 187 DO l=1,llm 188 zvert(l)=1. 189 ENDDO 190 fact=2. 191 DO l = 1, llm 192 zz = 1. - preff/presnivs(l) 193 zvert(l)= fact -( fact-1.)/( 1.+zz*zz ) 194 ENDDO 195 endif 196 197 198 PRINT*,'Constantes de temps de la diffusion horizontale' 199 200 tetamin = 1.e+6 201 202 DO l=1,llm 203 tetaudiv(l) = zvert(l)/tetagdiv 204 tetaurot(l) = zvert(l)/tetagrot 205 tetah(l) = zvert(l)/tetatemp 206 207 IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 208 IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 209 IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) 210 ENDDO 211 212 PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod 213 idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod 214 PRINT *,' INIDI tetamin idissip ',tetamin,idissip 215 idissip = MAX(iperiod,idissip) 216 dtdiss = idissip * dtvr 217 PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss 218 219 DO l = 1,llm 220 PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), 221 * dtdiss*tetah(l) 222 ENDDO 223 224 c 225 RETURN 226 END 53 ENDDO 54 ENDDO 55 56 idum = -1 57 zh(1) = RAN1(idum)-.5 58 idum = 0 59 DO ij = 2, ip1jmp1 60 zh(ij) = RAN1(idum) -.5 61 ENDDO 62 63 CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) 64 65 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 66 67 IF ( zhmin .GE. zhmax ) THEN 68 PRINT*,' Inidissip zh min max ',zhmin,zhmax 69 STOP'probleme generateur alleatoire dans inidissip' 70 ENDIF 71 72 zllm = ABS( zhmax ) 73 DO l = 1,50 74 IF(lstardis) THEN 75 CALL divgrad2(1,zh,deltap,niterh,zh) 76 ELSE 77 CALL divgrad (1,zh,niterh,zh) 78 ENDIF 79 80 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 81 82 zllm = ABS( zhmax ) 83 z1llm = 1./zllm 84 DO ij = 1,ip1jmp1 85 zh(ij) = zh(ij)* z1llm 86 ENDDO 87 ENDDO 88 89 IF(lstardis) THEN 90 cdivh = 1./ zllm 91 ELSE 92 cdivh = zllm ** ( -1./niterh ) 93 ENDIF 94 95 ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) 96 ! ----------------------------------------------------------------- 97 print*,'calcul des valeurs propres' 98 99 DO ii = 1, 2 100 ! 101 DO ij = 1, ip1jmp1 102 zu(ij) = RAN1(idum) -.5 103 ENDDO 104 CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) 105 DO ij = 1, ip1jm 106 zv(ij) = RAN1(idum) -.5 107 ENDDO 108 CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) 109 110 CALL minmax(iip1*jjp1,zu,umin,ullm ) 111 CALL minmax(iip1*jjm, zv,vmin,vllm ) 112 113 ullm = ABS ( ullm ) 114 vllm = ABS ( vllm ) 115 116 DO l = 1, 50 117 IF(ii.EQ.1) THEN 118 !cccc CALL covcont( 1,zu,zv,zu,zv ) 119 IF(lstardis) THEN 120 CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv ) 121 ELSE 122 CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv ) 123 ENDIF 124 ELSE 125 IF(lstardis) THEN 126 CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv ) 127 ELSE 128 CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv ) 129 ENDIF 130 ENDIF 131 132 CALL minmax(iip1*jjp1,zu,umin,ullm ) 133 CALL minmax(iip1*jjm, zv,vmin,vllm ) 134 135 ullm = ABS ( ullm ) 136 vllm = ABS ( vllm ) 137 138 zllm = MAX( ullm,vllm ) 139 z1llm = 1./ zllm 140 DO ij = 1, ip1jmp1 141 zu(ij) = zu(ij)* z1llm 142 ENDDO 143 DO ij = 1, ip1jm 144 zv(ij) = zv(ij)* z1llm 145 ENDDO 146 end DO 147 148 IF ( ii.EQ.1 ) THEN 149 IF(lstardis) THEN 150 cdivu = 1./zllm 151 ELSE 152 cdivu = zllm **( -1./nitergdiv ) 153 ENDIF 154 ELSE 155 IF(lstardis) THEN 156 crot = 1./ zllm 157 ELSE 158 crot = zllm **( -1./nitergrot ) 159 ENDIF 160 ENDIF 161 162 end DO 163 164 ! petit test pour les operateurs non star: 165 ! ---------------------------------------- 166 167 ! IF(.NOT.lstardis) THEN 168 fact = rad*24./REAL(jjm) 169 fact = fact*fact 170 PRINT*,'coef u ', fact/cdivu, 1./cdivu 171 PRINT*,'coef r ', fact/crot , 1./crot 172 PRINT*,'coef h ', fact/cdivh, 1./cdivh 173 ! ENDIF 174 175 !----------------------------------------------------------------------- 176 ! variation verticale du coefficient de dissipation: 177 ! -------------------------------------------------- 178 179 if (ok_strato .and. llm==39) then 180 do l=1,llm 181 pseudoz=8.*log(preff/presnivs(l)) 182 zvert(l)=1+ & 183 (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. & 184 *(dissip_factz-1.) 185 enddo 186 else 187 DO l=1,llm 188 zvert(l)=1. 189 ENDDO 190 fact=2. 191 DO l = 1, llm 192 zz = 1. - preff/presnivs(l) 193 zvert(l)= fact -( fact-1.)/( 1.+zz*zz ) 194 ENDDO 195 endif 196 197 198 PRINT*,'Constantes de temps de la diffusion horizontale' 199 200 tetamin = 1.e+6 201 202 DO l=1,llm 203 tetaudiv(l) = zvert(l)/tetagdiv 204 tetaurot(l) = zvert(l)/tetagrot 205 tetah(l) = zvert(l)/tetatemp 206 207 IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 208 IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 209 IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) 210 ENDDO 211 212 PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod 213 idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod 214 PRINT *,' INIDI tetamin idissip ',tetamin,idissip 215 idissip = MAX(iperiod,idissip) 216 dtdiss = idissip * dtvr 217 PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss 218 219 DO l = 1,llm 220 PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), & 221 dtdiss*tetah(l) 222 ENDDO 223 224 END SUBROUTINE inidissip
Note: See TracChangeset
for help on using the changeset viewer.
