- Timestamp:
- Jul 24, 2024, 2:54:37 PM (2 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inidissip.F90
r5106 r5116 1 2 1 ! $Id$ 3 2 4 SUBROUTINE inidissip( lstardis,nitergdiv,nitergrot,niterh, &5 tetagdiv,tetagrot,tetatemp, vert_prof_dissip)3 SUBROUTINE inidissip(lstardis, nitergdiv, nitergrot, niterh, & 4 tetagdiv, tetagrot, tetatemp, vert_prof_dissip) 6 5 !======================================================================= 7 6 ! initialisation de la dissipation horizontale … … 11 10 ! ------------- 12 11 13 USE control_mod, ONLY: dissip_period, iperiod12 USE control_mod, ONLY: dissip_period, iperiod 14 13 USE comconst_mod, ONLY: dissip_deltaz, dissip_factz, dissip_zref, & 15 14 dtdiss, dtvr, rad 16 15 USE comvert_mod, ONLY: preff, presnivs 17 16 USE lmdz_filtreg, ONLY: filtreg 17 USE lmdz_libmath, ONLY: minmax 18 18 19 19 IMPLICIT NONE … … 23 23 include "iniprint.h" 24 24 25 LOGICAL, INTENT(in) :: lstardis26 INTEGER, INTENT(in) :: nitergdiv,nitergrot,niterh27 REAL, INTENT(in) :: tetagdiv,tetagrot,tetatemp28 29 integer, INTENT(in) :: vert_prof_dissip25 LOGICAL, INTENT(in) :: lstardis 26 INTEGER, INTENT(in) :: nitergdiv, nitergrot, niterh 27 REAL, INTENT(in) :: tetagdiv, tetagrot, tetatemp 28 29 integer, INTENT(in) :: vert_prof_dissip 30 30 ! Vertical profile of horizontal dissipation 31 31 ! Allowed values: … … 33 33 ! 1: tanh of altitude 34 34 35 ! Local variables:36 REAL fact, zvert(llm),zz37 REAL zh(ip1jmp1), zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)38 real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1, llm)39 REAL ullm, vllm,umin,vmin,zhmin,zhmax35 ! Local variables: 36 REAL fact, zvert(llm), zz 37 REAL zh(ip1jmp1), zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1) 38 real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1, llm) 39 REAL ullm, vllm, umin, vmin, zhmin, zhmax 40 40 REAL zllm 41 41 42 INTEGER l, ij,idum,ii42 INTEGER l, ij, idum, ii 43 43 REAL tetamin 44 44 REAL pseudoz 45 character (len =80) :: abort_message45 character (len = 80) :: abort_message 46 46 47 47 REAL ran1 … … 53 53 ! ----------------------------------------------------------------- 54 54 55 crot 56 cdivu 57 cdivh 55 crot = -1. 56 cdivu = -1. 57 cdivh = -1. 58 58 59 59 ! calcul de la valeur propre de divgrad: … … 61 61 idum = 0 62 62 DO l = 1, llm 63 64 deltap(ij,l) = 1.65 66 ENDDO 67 68 idum 69 zh(1) = RAN1(idum) -.570 idum 63 DO ij = 1, ip1jmp1 64 deltap(ij, l) = 1. 65 ENDDO 66 ENDDO 67 68 idum = -1 69 zh(1) = RAN1(idum) - .5 70 idum = 0 71 71 DO ij = 2, ip1jmp1 72 zh(ij) = RAN1(idum) -.573 ENDDO 74 75 CALL filtreg (zh, jjp1,1,2,1,.TRUE.,1)76 77 CALL minmax(iip1 *jjp1,zh,zhmin,zhmax)78 79 IF ( zhmin >= zhmax) THEN80 write(lunout,*)' Inidissip zh min max ',zhmin,zhmax81 abort_message='probleme generateur alleatoire dans inidissip'82 CALL abort_gcm('inidissip',abort_message,1)72 zh(ij) = RAN1(idum) - .5 73 ENDDO 74 75 CALL filtreg (zh, jjp1, 1, 2, 1, .TRUE., 1) 76 77 CALL minmax(iip1 * jjp1, zh, zhmin, zhmax) 78 79 IF (zhmin >= zhmax) THEN 80 WRITE(lunout, *)' Inidissip zh min max ', zhmin, zhmax 81 abort_message = 'probleme generateur alleatoire dans inidissip' 82 CALL abort_gcm('inidissip', abort_message, 1) 83 83 ENDIF 84 84 85 zllm = ABS( zhmax)86 DO l = 1, 5087 88 CALL divgrad2(1,zh,deltap,niterh,divgra)89 90 CALL divgrad (1,zh,niterh,divgra)91 92 93 zllm= ABS(maxval(divgra))94 85 zllm = ABS(zhmax) 86 DO l = 1, 50 87 IF(lstardis) THEN 88 CALL divgrad2(1, zh, deltap, niterh, divgra) 89 ELSE 90 CALL divgrad (1, zh, niterh, divgra) 91 ENDIF 92 93 zllm = ABS(maxval(divgra)) 94 zh = divgra / zllm 95 95 ENDDO 96 96 97 97 IF(lstardis) THEN 98 cdivh = 1./ zllm98 cdivh = 1. / zllm 99 99 ELSE 100 cdivh = zllm ** ( -1./niterh)100 cdivh = zllm ** (-1. / niterh) 101 101 ENDIF 102 102 103 103 ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) 104 104 ! ----------------------------------------------------------------- 105 write(lunout,*)'inidissip: calcul des valeurs propres'105 WRITE(lunout, *)'inidissip: calcul des valeurs propres' 106 106 107 107 DO ii = 1, 2 108 108 109 DO ij = 1, ip1jmp1 110 zu(ij) = RAN1(idum) -.5 111 ENDDO 112 CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) 113 DO ij = 1, ip1jm 114 zv(ij) = RAN1(idum) -.5 115 ENDDO 116 CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) 117 118 CALL minmax(iip1*jjp1,zu,umin,ullm ) 119 CALL minmax(iip1*jjm, zv,vmin,vllm ) 120 121 ullm = ABS ( ullm ) 122 vllm = ABS ( vllm ) 123 124 DO l = 1, 50 125 IF(ii==1) THEN 126 !cccc CALL covcont( 1,zu,zv,zu,zv ) 127 IF(lstardis) THEN 128 CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy ) 129 ELSE 130 CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy ) 131 ENDIF 109 DO ij = 1, ip1jmp1 110 zu(ij) = RAN1(idum) - .5 111 ENDDO 112 CALL filtreg (zu, jjp1, 1, 2, 1, .TRUE., 1) 113 DO ij = 1, ip1jm 114 zv(ij) = RAN1(idum) - .5 115 ENDDO 116 CALL filtreg (zv, jjm, 1, 2, 1, .FALSE., 1) 117 118 CALL minmax(iip1 * jjp1, zu, umin, ullm) 119 CALL minmax(iip1 * jjm, zv, vmin, vllm) 120 121 ullm = ABS (ullm) 122 vllm = ABS (vllm) 123 124 DO l = 1, 50 125 IF(ii==1) THEN 126 !cccc CALL covcont( 1,zu,zv,zu,zv ) 127 IF(lstardis) THEN 128 CALL gradiv2(1, zu, zv, nitergdiv, gx, gy) 132 129 ELSE 133 IF(lstardis) THEN 134 CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy ) 135 ELSE 136 CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy ) 137 ENDIF 130 CALL gradiv (1, zu, zv, nitergdiv, gx, gy) 138 131 ENDIF 139 140 zllm = max(abs(maxval(gx)), abs(maxval(gy))) 141 zu = gx / zllm 142 zv = gy / zllm 143 end DO 144 145 IF ( ii==1 ) THEN 132 ELSE 146 133 IF(lstardis) THEN 147 cdivu = 1./zllm134 CALL nxgraro2(1, zu, zv, nitergrot, gx, gy) 148 135 ELSE 149 cdivu = zllm **( -1./nitergdiv)136 CALL nxgrarot(1, zu, zv, nitergrot, gx, gy) 150 137 ENDIF 151 ELSE 152 IF(lstardis) THEN 153 crot = 1./ zllm 154 ELSE 155 crot = zllm **( -1./nitergrot ) 156 ENDIF 157 ENDIF 138 ENDIF 139 140 zllm = max(abs(maxval(gx)), abs(maxval(gy))) 141 zu = gx / zllm 142 zv = gy / zllm 143 end DO 144 145 IF (ii==1) THEN 146 IF(lstardis) THEN 147 cdivu = 1. / zllm 148 ELSE 149 cdivu = zllm **(-1. / nitergdiv) 150 ENDIF 151 ELSE 152 IF(lstardis) THEN 153 crot = 1. / zllm 154 ELSE 155 crot = zllm **(-1. / nitergrot) 156 ENDIF 157 ENDIF 158 158 159 159 end DO … … 163 163 164 164 ! IF(.NOT.lstardis) THEN 165 fact = rad*24./REAL(jjm)166 fact = fact*fact167 write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu168 write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot169 write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh165 fact = rad * 24. / REAL(jjm) 166 fact = fact * fact 167 WRITE(lunout, *)'inidissip: coef u ', fact / cdivu, 1. / cdivu 168 WRITE(lunout, *)'inidissip: coef r ', fact / crot, 1. / crot 169 WRITE(lunout, *)'inidissip: coef h ', fact / cdivh, 1. / cdivh 170 170 ! ENDIF 171 171 … … 174 174 ! -------------------------------------------------- 175 175 176 if (vert_prof_dissip == 1) then177 do l=1,llm178 pseudoz=8.*log(preff/presnivs(l))179 zvert(l)=1+ &180 (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &181 *(dissip_factz-1.)182 176 if (vert_prof_dissip == 1) THEN 177 do l = 1, llm 178 pseudoz = 8. * log(preff / presnivs(l)) 179 zvert(l) = 1 + & 180 (tanh((pseudoz - dissip_zref) / dissip_deltaz) + 1.) / 2. & 181 * (dissip_factz - 1.) 182 enddo 183 183 else 184 DO l=1,llm185 zvert(l)=1.186 187 fact=2.188 189 zz = 1. - preff/presnivs(l)190 zvert(l)= fact -( fact-1.)/( 1.+zz*zz)191 184 DO l = 1, llm 185 zvert(l) = 1. 186 ENDDO 187 fact = 2. 188 DO l = 1, llm 189 zz = 1. - preff / presnivs(l) 190 zvert(l) = fact - (fact - 1.) / (1. + zz * zz) 191 ENDDO 192 192 endif 193 193 194 195 write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' 196 197 tetamin = 1.e+6 198 199 DO l=1,llm 200 tetaudiv(l) = zvert(l)/tetagdiv 201 tetaurot(l) = zvert(l)/tetagrot 202 tetah(l) = zvert(l)/tetatemp 203 204 IF( tetamin> (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 205 IF( tetamin> (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 206 IF( tetamin> (1./ tetah(l)) ) tetamin = 1./ tetah(l) 194 WRITE(lunout, *)'inidissip: Constantes de temps de la diffusion horizontale' 195 196 tetamin = 1.e+6 197 198 DO l = 1, llm 199 tetaudiv(l) = zvert(l) / tetagdiv 200 tetaurot(l) = zvert(l) / tetagrot 201 tetah(l) = zvert(l) / tetatemp 202 203 IF(tetamin> (1. / tetaudiv(l))) tetamin = 1. / tetaudiv(l) 204 IF(tetamin> (1. / tetaurot(l))) tetamin = 1. / tetaurot(l) 205 IF(tetamin> (1. / tetah(l))) tetamin = 1. / tetah(l) 207 206 ENDDO 208 207 209 208 ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def 210 209 IF (dissip_period == 0) THEN 211 dissip_period = INT( tetamin/( 2.*dtvr*iperiod)) * iperiod212 write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period213 dissip_period = MAX(iperiod,dissip_period)210 dissip_period = INT(tetamin / (2. * dtvr * iperiod)) * iperiod 211 WRITE(lunout, *)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ', tetamin, dtvr, iperiod, dissip_period 212 dissip_period = MAX(iperiod, dissip_period) 214 213 END IF 215 214 216 dtdiss 217 write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr218 219 DO l = 1, llm220 write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &221 dtdiss*tetah(l)215 dtdiss = dissip_period * dtvr 216 WRITE(lunout, *)'inidissip: dissip_period=', dissip_period, ' dtdiss=', dtdiss, ' dtvr=', dtvr 217 218 DO l = 1, llm 219 WRITE(lunout, *)zvert(l), dtdiss * tetaudiv(l), dtdiss * tetaurot(l), & 220 dtdiss * tetah(l) 222 221 ENDDO 223 222
Note: See TracChangeset
for help on using the changeset viewer.
