Changeset 270 for trunk/LMDZ.COMMON/libf/dyn3d
- Timestamp:
- Aug 18, 2011, 12:09:27 PM (13 years ago)
- Location:
- trunk/LMDZ.COMMON/libf/dyn3d
- Files:
-
- 1 added
- 14 edited
- 2 moved
-
advtrac.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F) (1 diff)
-
ce0l.F90 (modified) (5 diffs)
-
comdissipn.h (modified) (1 diff)
-
conf_gcm.F (modified) (5 diffs)
-
control_mod.F90 (modified) (1 diff)
-
defrun.F (modified) (1 diff)
-
gcm.F (modified) (2 diffs)
-
grilles_gcm_netcdf_sub.F90 (added)
-
ini_paramLMDZ_dyn.h (modified) (1 diff)
-
iniacademic.F90 (modified) (3 diffs)
-
iniconst.F (modified) (1 diff)
-
inidissip.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F) (1 diff)
-
integrd.F (modified) (2 diffs)
-
leapfrog.F (modified) (5 diffs)
-
limit_netcdf.F90 (modified) (22 diffs)
-
vlsplt.F (modified) (3 diffs)
-
write_paramLMDZ_dyn.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90
r269 r270 1 ! 2 ! $Id: advtrac.F 1446 2010-10-22 09:27:25Z emillour $ 3 ! 4 c 5 c 6 SUBROUTINE advtrac(pbaru,pbarv , 7 * p, masse,q,iapptrac,teta, 8 * flxw, 9 * pk) 10 c Auteur : F. Hourdin 11 c 12 c Modif. P. Le Van (20/12/97) 13 c F. Codron (10/99) 14 c D. Le Croller (07/2001) 15 c M.A Filiberti (04/2002) 16 c 17 USE infotrac 18 USE control_mod 19 20 21 IMPLICIT NONE 22 c 23 #include "dimensions.h" 24 #include "paramet.h" 25 #include "comconst.h" 26 #include "comvert.h" 27 #include "comdissip.h" 28 #include "comgeom2.h" 29 #include "logic.h" 30 #include "temps.h" 31 #include "ener.h" 32 #include "description.h" 33 #include "iniprint.h" 34 35 c------------------------------------------------------------------- 36 c Arguments 37 c------------------------------------------------------------------- 38 c Ajout PPM 39 c-------------------------------------------------------- 40 REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm) 41 c-------------------------------------------------------- 42 INTEGER iapptrac 43 REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) 44 REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm) 45 REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm) 46 REAL pk(ip1jmp1,llm) 47 REAL flxw(ip1jmp1,llm) 48 49 c------------------------------------------------------------- 50 c Variables locales 51 c------------------------------------------------------------- 52 53 REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) 54 REAL massem(ip1jmp1,llm),zdp(ip1jmp1) 55 REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 56 REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu 57 INTEGER iadvtr 58 INTEGER ij,l,iq,iiq 59 REAL zdpmin, zdpmax 60 EXTERNAL minmax 61 SAVE iadvtr, massem, pbaruc, pbarvc 62 DATA iadvtr/0/ 63 c---------------------------------------------------------- 64 c Rajouts pour PPM 65 c---------------------------------------------------------- 66 INTEGER indice,n 67 REAL dtbon ! Pas de temps adaptatif pour que CFL<1 68 REAL CFLmaxz,aaa,bbb ! CFL maximum 69 REAL psppm(iim,jjp1) ! pression au sol 70 REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm) 71 REAL qppm(iim*jjp1,llm,nqtot) 72 REAL fluxwppm(iim,jjp1,llm) 73 REAL apppm(llmp1), bpppm(llmp1) 74 LOGICAL dum,fill 75 DATA fill/.true./ 76 DATA dum/.true./ 77 78 integer,save :: countcfl=0 79 real cflx(ip1jmp1,llm) 80 real cfly(ip1jm,llm) 81 real cflz(ip1jmp1,llm) 82 real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm) 83 84 IF(iadvtr.EQ.0) THEN 85 CALL initial0(ijp1llm,pbaruc) 86 CALL initial0(ijmllm,pbarvc) 87 ENDIF 88 89 c accumulation des flux de masse horizontaux 90 DO l=1,llm 91 DO ij = 1,ip1jmp1 92 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 93 ENDDO 94 DO ij = 1,ip1jm 95 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 96 ENDDO 97 ENDDO 98 99 c selection de la masse instantannee des mailles avant le transport. 100 IF(iadvtr.EQ.0) THEN 101 102 CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) 103 ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 104 c 105 ENDIF 106 107 iadvtr = iadvtr+1 108 iapptrac = iadvtr 109 110 111 c Test pour savoir si on advecte a ce pas de temps 112 IF ( iadvtr.EQ.iapp_tracvl ) THEN 113 114 cc .. Modif P.Le Van ( 20/12/97 ) .... 115 cc 116 117 c traitement des flux de masse avant advection. 118 c 1. calcul de w 119 c 2. groupement des mailles pres du pole. 120 121 CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) 122 123 ! ... Flux de masse diaganostiques traceurs 124 flxw = wg / REAL(iapp_tracvl) 125 126 c test sur l'eventuelle creation de valeurs negatives de la masse 127 DO l=1,llm-1 128 DO ij = iip2+1,ip1jm 129 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) 130 s - pbarvg(ij-iip1,l) + pbarvg(ij,l) 131 s + wg(ij,l+1) - wg(ij,l) 132 ENDDO 133 CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) 134 DO ij = iip2,ip1jm 135 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 136 ENDDO 137 138 139 CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) 140 141 IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN 142 PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, 143 s ' MAX:', zdpmax 144 ENDIF 145 146 ENDDO 147 148 149 c------------------------------------------------------------------- 150 ! Calcul des criteres CFL en X, Y et Z 151 c------------------------------------------------------------------- 152 153 if (countcfl == 0. ) then 154 cflxmax(:)=0. 155 cflymax(:)=0. 156 cflzmax(:)=0. 157 endif 158 159 countcfl=countcfl+iapp_tracvl 160 cflx(:,:)=0. 161 cfly(:,:)=0. 162 cflz(:,:)=0. 163 do l=1,llm 164 do ij=iip2,ip1jm-1 165 if (pbarug(ij,l)>=0.) then 166 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) 167 else 168 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) 169 endif 170 enddo 171 enddo 172 do l=1,llm 173 do ij=iip2,ip1jm-1,iip1 174 cflx(ij+iip1,l)=cflx(ij,l) 175 enddo 176 enddo 177 178 do l=1,llm 179 do ij=1,ip1jm 180 if (pbarvg(ij,l)>=0.) then 181 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) 182 else 183 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) 184 endif 185 enddo 186 enddo 187 188 do l=2,llm 189 do ij=1,ip1jm 190 if (wg(ij,l)>=0.) then 191 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) 192 else 193 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) 194 endif 195 enddo 196 enddo 197 198 do l=1,llm 199 cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) 200 cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) 201 cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) 202 enddo 1 ! $Id: advtrac.F90 1549 2011-07-05 08:41:12Z lguez $ 2 3 SUBROUTINE advtrac(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk) 4 ! Auteur : F. Hourdin 5 ! 6 ! Modif. P. Le Van (20/12/97) 7 ! F. Codron (10/99) 8 ! D. Le Croller (07/2001) 9 ! M.A Filiberti (04/2002) 10 ! 11 USE infotrac 12 USE control_mod 13 14 15 IMPLICIT NONE 16 ! 17 include "dimensions.h" 18 include "paramet.h" 19 include "comconst.h" 20 include "comvert.h" 21 include "comdissip.h" 22 include "comgeom2.h" 23 include "logic.h" 24 include "temps.h" 25 include "ener.h" 26 include "description.h" 27 include "iniprint.h" 28 29 !------------------------------------------------------------------- 30 ! Arguments 31 !------------------------------------------------------------------- 32 ! Ajout PPM 33 !-------------------------------------------------------- 34 REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm) 35 !-------------------------------------------------------- 36 INTEGER iapptrac 37 REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) 38 REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm) 39 REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm) 40 REAL pk(ip1jmp1,llm) 41 REAL flxw(ip1jmp1,llm) 42 43 !------------------------------------------------------------- 44 ! Variables locales 45 !------------------------------------------------------------- 46 47 REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) 48 REAL massem(ip1jmp1,llm),zdp(ip1jmp1) 49 REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 50 REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu 51 INTEGER iadvtr 52 INTEGER ij,l,iq,iiq 53 REAL zdpmin, zdpmax 54 EXTERNAL minmax 55 SAVE iadvtr, massem, pbaruc, pbarvc 56 DATA iadvtr/0/ 57 !---------------------------------------------------------- 58 ! Rajouts pour PPM 59 !---------------------------------------------------------- 60 INTEGER indice,n 61 REAL dtbon ! Pas de temps adaptatif pour que CFL<1 62 REAL CFLmaxz,aaa,bbb ! CFL maximum 63 REAL psppm(iim,jjp1) ! pression au sol 64 REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm) 65 REAL qppm(iim*jjp1,llm,nqtot) 66 REAL fluxwppm(iim,jjp1,llm) 67 REAL apppm(llmp1), bpppm(llmp1) 68 LOGICAL dum,fill 69 DATA fill/.true./ 70 DATA dum/.true./ 71 72 integer,save :: countcfl=0 73 real cflx(ip1jmp1,llm) 74 real cfly(ip1jm,llm) 75 real cflz(ip1jmp1,llm) 76 real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm) 77 78 IF(iadvtr.EQ.0) THEN 79 CALL initial0(ijp1llm,pbaruc) 80 CALL initial0(ijmllm,pbarvc) 81 ENDIF 82 83 ! accumulation des flux de masse horizontaux 84 DO l=1,llm 85 DO ij = 1,ip1jmp1 86 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 87 ENDDO 88 DO ij = 1,ip1jm 89 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 90 ENDDO 91 ENDDO 92 93 ! selection de la masse instantannee des mailles avant le transport. 94 IF(iadvtr.EQ.0) THEN 95 96 CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) 97 !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 98 ! 99 ENDIF 100 101 iadvtr = iadvtr+1 102 iapptrac = iadvtr 103 104 105 ! Test pour savoir si on advecte a ce pas de temps 106 IF ( iadvtr.EQ.iapp_tracvl ) THEN 107 108 !c .. Modif P.Le Van ( 20/12/97 ) .... 109 !c 110 111 ! traitement des flux de masse avant advection. 112 ! 1. calcul de w 113 ! 2. groupement des mailles pres du pole. 114 115 CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) 116 117 ! ... Flux de masse diaganostiques traceurs 118 flxw = wg / REAL(iapp_tracvl) 119 120 ! test sur l'eventuelle creation de valeurs negatives de la masse 121 DO l=1,llm-1 122 DO ij = iip2+1,ip1jm 123 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) & 124 - pbarvg(ij-iip1,l) + pbarvg(ij,l) & 125 + wg(ij,l+1) - wg(ij,l) 126 ENDDO 127 CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) 128 DO ij = iip2,ip1jm 129 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 130 ENDDO 131 132 133 CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) 134 135 IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN 136 PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, & 137 ' MAX:', zdpmax 138 ENDIF 139 140 ENDDO 141 142 143 !------------------------------------------------------------------- 144 ! Calcul des criteres CFL en X, Y et Z 145 !------------------------------------------------------------------- 146 147 if (countcfl == 0. ) then 148 cflxmax(:)=0. 149 cflymax(:)=0. 150 cflzmax(:)=0. 151 endif 152 153 countcfl=countcfl+iapp_tracvl 154 cflx(:,:)=0. 155 cfly(:,:)=0. 156 cflz(:,:)=0. 157 do l=1,llm 158 do ij=iip2,ip1jm-1 159 if (pbarug(ij,l)>=0.) then 160 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) 161 else 162 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) 163 endif 164 enddo 165 enddo 166 do l=1,llm 167 do ij=iip2,ip1jm-1,iip1 168 cflx(ij+iip1,l)=cflx(ij,l) 169 enddo 170 enddo 171 172 do l=1,llm 173 do ij=1,ip1jm 174 if (pbarvg(ij,l)>=0.) then 175 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) 176 else 177 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) 178 endif 179 enddo 180 enddo 181 182 do l=2,llm 183 do ij=1,ip1jm 184 if (wg(ij,l)>=0.) then 185 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) 186 else 187 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) 188 endif 189 enddo 190 enddo 191 192 do l=1,llm 193 cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) 194 cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) 195 cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) 196 enddo 203 197 204 198 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 205 ! Par defaut, on sort le diagnostic des CFL tous les jours.206 ! Si on veut le sortir a chaque pas d'advection en cas de plantage207 ! if (countcfl==iapp_tracvl) then199 ! Par defaut, on sort le diagnostic des CFL tous les jours. 200 ! Si on veut le sortir a chaque pas d'advection en cas de plantage 201 ! if (countcfl==iapp_tracvl) then 208 202 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 209 210 211 write(lunout,*) 'L, CFLmax '212 s,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))213 214 215 216 217 c-------------------------------------------------------------------218 cAdvection proprement dite (Modification Le Croller (07/2001)219 c-------------------------------------------------------------------220 221 c----------------------------------------------------222 cCalcul des moyennes basées sur la masse223 c----------------------------------------------------224 225 226 c-----------------------------------------------------------227 cAppel des sous programmes d'advection228 c-----------------------------------------------------------229 230 ccall clock(t_initial)203 if (countcfl==day_step) then 204 do l=1,llm 205 write(lunout,*) 'L, CFLmax ' & 206 ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l)) 207 enddo 208 countcfl=0 209 endif 210 211 !------------------------------------------------------------------- 212 ! Advection proprement dite (Modification Le Croller (07/2001) 213 !------------------------------------------------------------------- 214 215 !---------------------------------------------------- 216 ! Calcul des moyennes basées sur la masse 217 !---------------------------------------------------- 218 call massbar(massem,massebx,masseby) 219 220 !----------------------------------------------------------- 221 ! Appel des sous programmes d'advection 222 !----------------------------------------------------------- 223 do iq=1,nqtot 224 ! call clock(t_initial) 231 225 if(iadv(iq) == 0) cycle 232 c----------------------------------------------------------------233 cSchema de Van Leer I MUSCL234 c----------------------------------------------------------------226 ! ---------------------------------------------------------------- 227 ! Schema de Van Leer I MUSCL 228 ! ---------------------------------------------------------------- 235 229 if(iadv(iq).eq.10) THEN 236 237 238 c----------------------------------------------------------------239 cSchema "pseudo amont" + test sur humidite specifique240 Cpour la vapeur d'eau. F. Codron241 c----------------------------------------------------------------230 call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr) 231 232 ! ---------------------------------------------------------------- 233 ! Schema "pseudo amont" + test sur humidite specifique 234 ! pour la vapeur d'eau. F. Codron 235 ! ---------------------------------------------------------------- 242 236 else if(iadv(iq).eq.14) then 243 c 244 CALL vlspltqs( q(1,1,1), 2., massem, wg , 245 *pbarug,pbarvg,dtvr,p,pk,teta )246 c----------------------------------------------------------------247 cSchema de Frederic Hourdin248 c----------------------------------------------------------------237 ! 238 CALL vlspltqs( q(1,1,1), 2., massem, wg , & 239 pbarug,pbarvg,dtvr,p,pk,teta ) 240 ! ---------------------------------------------------------------- 241 ! Schema de Frederic Hourdin 242 ! ---------------------------------------------------------------- 249 243 else if(iadv(iq).eq.12) then 250 cPas de temps adaptatif244 ! Pas de temps adaptatif 251 245 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 252 246 if (n.GT.1) then 253 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',254 sdtvr,'n=',n247 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 248 dtvr,'n=',n 255 249 endif 256 250 do indice=1,n 257 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)251 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) 258 252 end do 259 253 else if(iadv(iq).eq.13) then 260 cPas de temps adaptatif254 ! Pas de temps adaptatif 261 255 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 262 256 if (n.GT.1) then 263 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',264 sdtvr,'n=',n265 endif 266 do indice=1,n267 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)268 end do269 c----------------------------------------------------------------270 cSchema de pente SLOPES271 c----------------------------------------------------------------257 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 258 dtvr,'n=',n 259 endif 260 do indice=1,n 261 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) 262 end do 263 ! ---------------------------------------------------------------- 264 ! Schema de pente SLOPES 265 ! ---------------------------------------------------------------- 272 266 else if (iadv(iq).eq.20) then 273 274 275 c----------------------------------------------------------------276 cSchema de Prather277 c----------------------------------------------------------------267 call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) 268 269 ! ---------------------------------------------------------------- 270 ! Schema de Prather 271 ! ---------------------------------------------------------------- 278 272 else if (iadv(iq).eq.30) then 279 cPas de temps adaptatif273 ! Pas de temps adaptatif 280 274 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 281 275 if (n.GT.1) then 282 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', 283 s dtvr,'n=',n 284 endif 285 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, 286 s n,dtbon) 287 288 c ---------------------------------------------------------------- 289 c Schemas PPM Lin et Rood 290 c ---------------------------------------------------------------- 291 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. 292 s iadv(iq).LE.18)) then 293 294 c Test sur le flux horizontal 295 c Pas de temps adaptatif 296 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 297 if (n.GT.1) then 298 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', 299 s dtvr,'n=',n 300 endif 301 c Test sur le flux vertical 302 CFLmaxz=0. 303 do l=2,llm 304 do ij=iip2,ip1jm 305 aaa=wg(ij,l)*dtvr/massem(ij,l) 306 CFLmaxz=max(CFLmaxz,aaa) 307 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 308 CFLmaxz=max(CFLmaxz,bbb) 276 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 277 dtvr,'n=',n 278 endif 279 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, & 280 n,dtbon) 281 282 ! ---------------------------------------------------------------- 283 ! Schemas PPM Lin et Rood 284 ! ---------------------------------------------------------------- 285 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. & 286 iadv(iq).LE.18)) then 287 288 ! Test sur le flux horizontal 289 ! Pas de temps adaptatif 290 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 291 if (n.GT.1) then 292 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 293 dtvr,'n=',n 294 endif 295 ! Test sur le flux vertical 296 CFLmaxz=0. 297 do l=2,llm 298 do ij=iip2,ip1jm 299 aaa=wg(ij,l)*dtvr/massem(ij,l) 300 CFLmaxz=max(CFLmaxz,aaa) 301 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 302 CFLmaxz=max(CFLmaxz,bbb) 303 enddo 309 304 enddo 310 enddo 311 if (CFLmaxz.GE.1) then 312 write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 313 endif 314 315 c----------------------------------------------------------- 316 c Ss-prg interface LMDZ.4->PPM3d 317 c----------------------------------------------------------- 318 319 call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, 320 s apppm,bpppm,massebx,masseby,pbarug,pbarvg, 321 s unatppm,vnatppm,psppm) 322 323 do indice=1,n 324 c--------------------------------------------------------------------- 325 c VL (version PPM) horiz. et PPM vert. 326 c--------------------------------------------------------------------- 327 if (iadv(iq).eq.11) then 328 c Ss-prg PPM3d de Lin 329 call ppm3d(1,qppm(1,1,iq), 330 s psppm,psppm, 331 s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, 332 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 333 s fill,dum,220.) 334 335 c---------------------------------------------------------------------- 336 c Monotonic PPM 337 c---------------------------------------------------------------------- 338 else if (iadv(iq).eq.16) then 339 c Ss-prg PPM3d de Lin 340 call ppm3d(1,qppm(1,1,iq), 341 s psppm,psppm, 342 s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, 343 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 344 s fill,dum,220.) 345 c--------------------------------------------------------------------- 346 347 c--------------------------------------------------------------------- 348 c Semi Monotonic PPM 349 c--------------------------------------------------------------------- 350 else if (iadv(iq).eq.17) then 351 c Ss-prg PPM3d de Lin 352 call ppm3d(1,qppm(1,1,iq), 353 s psppm,psppm, 354 s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, 355 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 356 s fill,dum,220.) 357 c--------------------------------------------------------------------- 358 359 c--------------------------------------------------------------------- 360 c Positive Definite PPM 361 c--------------------------------------------------------------------- 362 else if (iadv(iq).eq.18) then 363 c Ss-prg PPM3d de Lin 364 call ppm3d(1,qppm(1,1,iq), 365 s psppm,psppm, 366 s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, 367 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 368 s fill,dum,220.) 369 c--------------------------------------------------------------------- 370 endif 371 enddo 372 c----------------------------------------------------------------- 373 c Ss-prg interface PPM3d-LMDZ.4 374 c----------------------------------------------------------------- 375 call interpost(q(1,1,iq),qppm(1,1,iq)) 376 endif 377 c---------------------------------------------------------------------- 378 379 c----------------------------------------------------------------- 380 c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 381 c et Nord j=1 382 c----------------------------------------------------------------- 383 384 c call traceurpole(q(1,1,iq),massem) 385 386 c calcul du temps cpu pour un schema donne 387 388 c call clock(t_final) 389 cym tps_cpu=t_final-t_initial 390 cym cpuadv(iq)=cpuadv(iq)+tps_cpu 391 392 end DO 393 394 395 c------------------------------------------------------------------ 396 c on reinitialise a zero les flux de masse cumules 397 c--------------------------------------------------- 398 iadvtr=0 399 400 ENDIF ! if iadvtr.EQ.iapp_tracvl 401 402 RETURN 403 END 404 305 if (CFLmaxz.GE.1) then 306 write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 307 endif 308 309 !----------------------------------------------------------- 310 ! Ss-prg interface LMDZ.4->PPM3d 311 !----------------------------------------------------------- 312 313 call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, & 314 apppm,bpppm,massebx,masseby,pbarug,pbarvg, & 315 unatppm,vnatppm,psppm) 316 317 do indice=1,n 318 !---------------------------------------------------------------- 319 ! VL (version PPM) horiz. et PPM vert. 320 !---------------------------------------------------------------- 321 if (iadv(iq).eq.11) then 322 ! Ss-prg PPM3d de Lin 323 call ppm3d(1,qppm(1,1,iq), & 324 psppm,psppm, & 325 unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, & 326 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 327 fill,dum,220.) 328 329 !------------------------------------------------------------- 330 ! Monotonic PPM 331 !------------------------------------------------------------- 332 else if (iadv(iq).eq.16) then 333 ! Ss-prg PPM3d de Lin 334 call ppm3d(1,qppm(1,1,iq), & 335 psppm,psppm, & 336 unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, & 337 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 338 fill,dum,220.) 339 !------------------------------------------------------------- 340 341 !------------------------------------------------------------- 342 ! Semi Monotonic PPM 343 !------------------------------------------------------------- 344 else if (iadv(iq).eq.17) then 345 ! Ss-prg PPM3d de Lin 346 call ppm3d(1,qppm(1,1,iq), & 347 psppm,psppm, & 348 unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, & 349 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 350 fill,dum,220.) 351 !------------------------------------------------------------- 352 353 !------------------------------------------------------------- 354 ! Positive Definite PPM 355 !------------------------------------------------------------- 356 else if (iadv(iq).eq.18) then 357 ! Ss-prg PPM3d de Lin 358 call ppm3d(1,qppm(1,1,iq), & 359 psppm,psppm, & 360 unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, & 361 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 362 fill,dum,220.) 363 !------------------------------------------------------------- 364 endif 365 enddo 366 !----------------------------------------------------------------- 367 ! Ss-prg interface PPM3d-LMDZ.4 368 !----------------------------------------------------------------- 369 call interpost(q(1,1,iq),qppm(1,1,iq)) 370 endif 371 !---------------------------------------------------------------------- 372 373 !----------------------------------------------------------------- 374 ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 375 ! et Nord j=1 376 !----------------------------------------------------------------- 377 378 ! call traceurpole(q(1,1,iq),massem) 379 380 ! calcul du temps cpu pour un schema donne 381 382 ! call clock(t_final) 383 !ym tps_cpu=t_final-t_initial 384 !ym cpuadv(iq)=cpuadv(iq)+tps_cpu 385 386 end DO 387 388 389 !------------------------------------------------------------------ 390 ! on reinitialise a zero les flux de masse cumules 391 !--------------------------------------------------- 392 iadvtr=0 393 394 ENDIF ! if iadvtr.EQ.iapp_tracvl 395 396 END SUBROUTINE advtrac -
trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90
r101 r270 1 1 ! 2 ! $Id: ce0l.F90 1 425 2010-09-02 13:45:23Z lguez$2 ! $Id: ce0l.F90 1511 2011-04-28 15:21:47Z jghattas $ 3 3 ! 4 4 !------------------------------------------------------------------------------- … … 27 27 #endif 28 28 IMPLICIT NONE 29 #include "iniprint.h"30 29 #ifndef CPP_EARTH 31 30 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' … … 37 36 #include "paramet.h" 38 37 #include "indicesol.h" 38 #include "iniprint.h" 39 39 #include "temps.h" 40 40 #include "logic.h" 41 INTEGER, PARAMETER :: longcles=20 42 REAL, DIMENSION(longcles) :: clesphy0 41 43 REAL, DIMENSION(iip1,jjp1) :: masque 42 44 CHARACTER(LEN=15) :: calnd 45 REAL, DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol 43 46 !------------------------------------------------------------------------------- 44 CALL conf_gcm( 99, .TRUE. )47 CALL conf_gcm( 99, .TRUE. , clesphy0 ) 45 48 46 49 CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) … … 78 81 WRITE(lunout,'(//)') 79 82 WRITE(lunout,*) ' interbar = ',interbar 80 CALL etat0_netcdf(interbar,masque, ok_etat0)83 CALL etat0_netcdf(interbar,masque,phis,ok_etat0) 81 84 82 85 IF(ok_limit) THEN … … 89 92 END IF 90 93 94 IF (grilles_gcm_netcdf) THEN 95 WRITE(lunout,'(//)') 96 WRITE(lunout,*) ' *************************** ' 97 WRITE(lunout,*) ' *** grilles_gcm_netcdf *** ' 98 WRITE(lunout,*) ' *************************** ' 99 WRITE(lunout,'(//)') 100 CALL grilles_gcm_netcdf_sub(masque,phis) 101 END IF 91 102 #endif 92 103 ! of #ifndef CPP_EARTH #else -
trunk/LMDZ.COMMON/libf/dyn3d/comdissipn.h
r1 r270 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 !----------------------------------------------------------------------- -
trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F
r119 r270 230 230 CALL getin('output_grads_dyn',output_grads_dyn) 231 231 232 !Config Key = idissip232 !Config Key = dissip_period 233 233 !Config Desc = periode de la dissipation 234 !Config Def = 10234 !Config Def = 0 235 235 !Config Help = periode de la dissipation 236 !Config (en pas) ... a completer ! 237 idissip = 10 238 CALL getin('idissip',idissip) 236 !Config dissip_period=0 => la valeur sera calcule dans inidissip 237 !Config dissip_period>0 => on prend cette valeur 238 dissip_period = 0 239 CALL getin('dissip_period',dissip_period) 239 240 240 241 ccc .... P. Le Van , modif le 29/04/97 .pour la dissipation ... … … 652 653 write(lunout,*)' periodav = ', periodav 653 654 write(lunout,*)' output_grads_dyn = ', output_grads_dyn 654 write(lunout,*)' idissip = ', idissip655 write(lunout,*)' dissip_period = ', dissip_period 655 656 write(lunout,*)' lstardis = ', lstardis 656 657 write(lunout,*)' nitergdiv = ', nitergdiv … … 876 877 CALL getin('ok_etat0',ok_etat0) 877 878 879 !Config Key = grilles_gcm_netcdf 880 !Config Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit 881 !Config Def = n 882 grilles_gcm_netcdf = .FALSE. 883 CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf) 884 878 885 write(lunout,*)' #########################################' 879 886 write(lunout,*)' Configuration des parametres de cel0' … … 890 897 write(lunout,*)' periodav = ', periodav 891 898 write(lunout,*)' output_grads_dyn = ', output_grads_dyn 892 write(lunout,*)' idissip = ', idissip899 write(lunout,*)' dissip_period = ', dissip_period 893 900 write(lunout,*)' lstardis = ', lstardis 894 901 write(lunout,*)' nitergdiv = ', nitergdiv … … 922 929 write(lunout,*)' ok_limit = ', ok_limit 923 930 write(lunout,*)' ok_etat0 = ', ok_etat0 931 write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf 924 932 c 925 933 RETURN -
trunk/LMDZ.COMMON/libf/dyn3d/control_mod.F90
r97 r270 12 12 REAL :: periodav 13 13 INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys 14 INTEGER :: iconser,iecri, idissip,iphysiq,iecrimoy14 INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy 15 15 INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn 16 16 LOGICAL :: offline -
trunk/LMDZ.COMMON/libf/dyn3d/defrun.F
r6 r270 139 139 140 140 READ (tapedef,9001) ch1,ch4 141 READ (tapedef,*) idissip142 WRITE(tapeout,9001) ch1,' idissip'143 WRITE(tapeout,*) idissip141 READ (tapedef,*) dissip_period 142 WRITE(tapeout,9001) ch1,'dissip_period' 143 WRITE(tapeout,*) dissip_period 144 144 145 145 ccc .... P. Le Van , modif le 29/04/97 .pour la dissipation ... -
trunk/LMDZ.COMMON/libf/dyn3d/gcm.F
r130 r270 427 427 c ------------------------------- 428 428 429 IF (call_iniphys.and.(iflag_phys .eq.1)) THEN429 IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN 430 430 latfi(1)=rlatu(1) 431 431 lonfi(1)=0. … … 486 486 #endif 487 487 488 #ifdef CPP_EARTH 489 ! Create start file (startphy.nc) and boundary conditions (limit.nc) 490 ! for the Earth verstion 491 if (iflag_phys>=100) then 492 call iniaqua(ngridmx,latfi,lonfi,iflag_phys) 493 endif 494 #endif 495 488 496 if (planet_type.eq."mars") then 489 497 ! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem0_mars -
trunk/LMDZ.COMMON/libf/dyn3d/ini_paramLMDZ_dyn.h
r1 r270 78 78 . "once", dt_cum,dt_cum) 79 79 c 80 CALL histdef(nid_ctesGCM, " idissip",80 CALL histdef(nid_ctesGCM, "dissip_period", 81 81 . "periode de la dissipation (en pas) ... a completer", 82 82 . "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, -
trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90
r127 r270 1 1 ! 2 ! $Id: iniacademic.F90 152 0 2011-05-23 11:37:09Z emillour$2 ! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $ 3 3 ! 4 4 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) … … 115 115 endif 116 116 117 academic_case: if (iflag_phys == 2) then117 academic_case: if (iflag_phys >= 2) then 118 118 ! initializations 119 119 … … 208 208 IF (.NOT. read_start) THEN 209 209 ! surface pressure 210 ps(:)=preff 210 if (iflag_phys>2) then 211 ps(:)=108080. ! Earth aqua/terra planets 212 else 213 ps(:)=preff 214 endif 211 215 ! ground geopotential 212 216 phis(:)=0. -
trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F
r127 r270 47 47 c----------------------------------------------------------------------- 48 48 49 dtdiss = idissip * dtvr50 49 dtphys = iphysiq * dtvr 51 50 unsim = 1./iim -
trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90
r269 r270 1 1 ! 2 ! $Id: inidissip.F 1403 2010-07-01 09:02:53Z fairhead$2 ! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $ 3 3 ! 4 SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , 5 * tetagdiv,tetagrot,tetatemp ) 6 c======================================================================= 7 c initialisation de la dissipation horizontale 8 c======================================================================= 9 c----------------------------------------------------------------------- 10 c declarations: 11 c ------------- 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 Pup 34 35 REAL ran1 36 37 38 c----------------------------------------------------------------------- 39 c 40 c calcul des valeurs propres des operateurs par methode iterrative: 41 c ----------------------------------------------------------------- 42 43 crot = -1. 44 cdivu = -1. 45 cdivh = -1. 46 47 c calcul de la valeur propre de divgrad: 48 c -------------------------------------- 49 idum = 0 50 DO l = 1, llm 51 DO ij = 1, ip1jmp1 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, only : dissip_period,iperiod 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 include "iniprint.h" 23 24 LOGICAL,INTENT(in) :: lstardis 25 INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh 26 REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp 27 28 ! Local variables: 29 REAL fact,zvert(llm),zz 30 REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm) 31 REAL ullm,vllm,umin,vmin,zhmin,zhmax 32 REAL zllm,z1llm 33 34 INTEGER l,ij,idum,ii 35 REAL tetamin 36 REAL Pup 37 character (len=80) :: abort_message 38 39 REAL ran1 40 41 42 !----------------------------------------------------------------------- 43 ! 44 ! calcul des valeurs propres des operateurs par methode iterrative: 45 ! ----------------------------------------------------------------- 46 47 crot = -1. 48 cdivu = -1. 49 cdivh = -1. 50 51 ! calcul de la valeur propre de divgrad: 52 ! -------------------------------------- 53 idum = 0 54 DO l = 1, llm 55 DO ij = 1, ip1jmp1 52 56 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 c First step: going from 1 to dissip_fac_mid (in gcm.def) 180 c============ 181 DO l=1,llm 182 zz = 1. - preff/presnivs(l) 183 zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) 184 ENDDO 185 186 write(*,*) 'Dissipation : ' 187 write(*,*) 'Multiplication de la dissipation en altitude :' 188 write(*,*) ' dissip_fac_mid =', dissip_fac_mid 189 190 c Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) 191 c========================== 192 c Utilisation de la fonction tangente hyperbolique pour augmenter 193 c arbitrairement la dissipation et donc la stabilite du modele a 194 c partir d'une certaine altitude. 195 196 c Le facteur multiplicatif de basse atmosphere etant deja pris 197 c en compte, il faut diviser le facteur multiplicatif de haute 198 c atmosphere par celui-ci. 199 200 if (ok_strato) then 201 202 Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) 203 do l=1,llm 204 zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) 205 & *(1-(0.5*(1+tanh(-6./dissip_deltaz* 206 & (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) 207 enddo 208 209 write(*,*) ' dissip_fac_up =', dissip_fac_up 210 write(*,*) 'Transition mid /up: Pupstart,delta =', 211 & dissip_pupstart,'Pa', dissip_deltaz , '(km)' 212 213 endif 214 215 216 PRINT*,'Constantes de temps de la diffusion horizontale' 217 218 tetamin = 1.e+6 219 220 DO l=1,llm 221 tetaudiv(l) = zvert(l)/tetagdiv 222 tetaurot(l) = zvert(l)/tetagrot 223 tetah(l) = zvert(l)/tetatemp 224 225 IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 226 IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 227 IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) 228 ENDDO 229 230 PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod 231 idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod 232 PRINT *,' INIDI tetamin idissip ',tetamin,idissip 233 idissip = MAX(iperiod,idissip) 234 dtdiss = idissip * dtvr 235 PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss 236 237 DO l = 1,llm 238 PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), 239 * dtdiss*tetah(l) 240 ENDDO 241 242 c 243 RETURN 244 END 57 ENDDO 58 ENDDO 59 60 idum = -1 61 zh(1) = RAN1(idum)-.5 62 idum = 0 63 DO ij = 2, ip1jmp1 64 zh(ij) = RAN1(idum) -.5 65 ENDDO 66 67 CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) 68 69 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 70 71 IF ( zhmin .GE. zhmax ) THEN 72 write(lunout,*)' Inidissip zh min max ',zhmin,zhmax 73 abort_message='probleme generateur alleatoire dans inidissip' 74 call abort_gcm('inidissip',abort_message,1) 75 ENDIF 76 77 zllm = ABS( zhmax ) 78 DO l = 1,50 79 IF(lstardis) THEN 80 CALL divgrad2(1,zh,deltap,niterh,zh) 81 ELSE 82 CALL divgrad (1,zh,niterh,zh) 83 ENDIF 84 85 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 86 87 zllm = ABS( zhmax ) 88 z1llm = 1./zllm 89 DO ij = 1,ip1jmp1 90 zh(ij) = zh(ij)* z1llm 91 ENDDO 92 ENDDO 93 94 IF(lstardis) THEN 95 cdivh = 1./ zllm 96 ELSE 97 cdivh = zllm ** ( -1./niterh ) 98 ENDIF 99 100 ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) 101 ! ----------------------------------------------------------------- 102 write(lunout,*)'inidissip: calcul des valeurs propres' 103 104 DO ii = 1, 2 105 ! 106 DO ij = 1, ip1jmp1 107 zu(ij) = RAN1(idum) -.5 108 ENDDO 109 CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) 110 DO ij = 1, ip1jm 111 zv(ij) = RAN1(idum) -.5 112 ENDDO 113 CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) 114 115 CALL minmax(iip1*jjp1,zu,umin,ullm ) 116 CALL minmax(iip1*jjm, zv,vmin,vllm ) 117 118 ullm = ABS ( ullm ) 119 vllm = ABS ( vllm ) 120 121 DO l = 1, 50 122 IF(ii.EQ.1) THEN 123 !cccc CALL covcont( 1,zu,zv,zu,zv ) 124 IF(lstardis) THEN 125 CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv ) 126 ELSE 127 CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv ) 128 ENDIF 129 ELSE 130 IF(lstardis) THEN 131 CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv ) 132 ELSE 133 CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv ) 134 ENDIF 135 ENDIF 136 137 CALL minmax(iip1*jjp1,zu,umin,ullm ) 138 CALL minmax(iip1*jjm, zv,vmin,vllm ) 139 140 ullm = ABS ( ullm ) 141 vllm = ABS ( vllm ) 142 143 zllm = MAX( ullm,vllm ) 144 z1llm = 1./ zllm 145 DO ij = 1, ip1jmp1 146 zu(ij) = zu(ij)* z1llm 147 ENDDO 148 DO ij = 1, ip1jm 149 zv(ij) = zv(ij)* z1llm 150 ENDDO 151 end DO 152 153 IF ( ii.EQ.1 ) THEN 154 IF(lstardis) THEN 155 cdivu = 1./zllm 156 ELSE 157 cdivu = zllm **( -1./nitergdiv ) 158 ENDIF 159 ELSE 160 IF(lstardis) THEN 161 crot = 1./ zllm 162 ELSE 163 crot = zllm **( -1./nitergrot ) 164 ENDIF 165 ENDIF 166 167 end DO 168 169 ! petit test pour les operateurs non star: 170 ! ---------------------------------------- 171 172 ! IF(.NOT.lstardis) THEN 173 fact = rad*24./REAL(jjm) 174 fact = fact*fact 175 write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu 176 write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot 177 write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh 178 ! ENDIF 179 180 !----------------------------------------------------------------------- 181 ! variation verticale du coefficient de dissipation: 182 ! -------------------------------------------------- 183 184 ! First step: going from 1 to dissip_fac_mid (in gcm.def) 185 !============ 186 DO l=1,llm 187 zz = 1. - preff/presnivs(l) 188 zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) 189 ENDDO 190 191 write(lunout,*) 'Dissipation : ' 192 write(lunout,*) 'Multiplication de la dissipation en altitude :' 193 write(lunout,*) ' dissip_fac_mid =', dissip_fac_mid 194 195 ! Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) 196 !========================== 197 ! Utilisation de la fonction tangente hyperbolique pour augmenter 198 ! arbitrairement la dissipation et donc la stabilite du modele a 199 ! partir d'une certaine altitude. 200 201 ! Le facteur multiplicatif de basse atmosphere etant deja pris 202 ! en compte, il faut diviser le facteur multiplicatif de haute 203 ! atmosphere par celui-ci. 204 205 if (ok_strato) then 206 207 Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) 208 do l=1,llm 209 zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) & 210 *(1-(0.5*(1+tanh(-6./dissip_deltaz* & 211 (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) 212 enddo 213 214 write(*,*) ' dissip_fac_up =', dissip_fac_up 215 write(*,*) 'Transition mid /up: Pupstart,delta =', & 216 dissip_pupstart,'Pa', dissip_deltaz , '(km)' 217 218 endif 219 220 221 write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' 222 223 tetamin = 1.e+6 224 225 DO l=1,llm 226 tetaudiv(l) = zvert(l)/tetagdiv 227 tetaurot(l) = zvert(l)/tetagrot 228 tetah(l) = zvert(l)/tetatemp 229 230 IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 231 IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 232 IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) 233 ENDDO 234 235 ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def 236 IF (dissip_period == 0) THEN 237 dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod 238 write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period 239 dissip_period = MAX(iperiod,dissip_period) 240 END IF 241 242 dtdiss = dissip_period * dtvr 243 write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr 244 245 DO l = 1,llm 246 write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), & 247 dtdiss*tetah(l) 248 ENDDO 249 250 END SUBROUTINE inidissip -
trunk/LMDZ.COMMON/libf/dyn3d/integrd.F
r7 r270 1 1 ! 2 ! $Id: integrd.F 1 446 2010-10-22 09:27:25Z emillour$2 ! $Id: integrd.F 1550 2011-07-05 09:44:55Z lguez $ 3 3 ! 4 4 SUBROUTINE integrd … … 89 89 IF( ps(ij).LT.0. ) THEN 90 90 PRINT*,' Au point ij = ',ij, ' , pression sol neg. ', ps(ij) 91 STOP' dans integrd' 91 print *, ' dans integrd' 92 stop 1 92 93 ENDIF 93 94 ENDDO -
trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F
r130 r270 152 152 logical ok_sync 153 153 parameter (ok_sync = .true.) 154 logical physic 154 155 155 156 data callinigrads/.true./ … … 222 223 223 224 itau = 0 225 physic=.true. 226 if (iflag_phys==0.or.iflag_phys==2) physic=.false. 227 224 228 c iday = day_ini+itau/day_step 225 229 c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0 … … 322 326 ! Purely Matsuno time stepping 323 327 IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE. 324 IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE. 328 IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward ) 329 s apdiss = .TRUE. 325 330 IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 326 s .and. iflag_phys.EQ.1) apphys = .TRUE.331 s .and. physic ) apphys = .TRUE. 327 332 ELSE 328 333 ! Leapfrog/Matsuno time stepping 329 334 IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE. 330 IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE. 331 IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE. 335 IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward ) 336 s apdiss = .TRUE. 337 IF( MOD(itau+1,iphysiq).EQ.0.AND.physic ) apphys=.TRUE. 332 338 END IF 333 339 … … 357 363 c ------------------------------------------------------------- 358 364 359 IF( forward. OR . leapf ) THEN 360 365 ! IF( forward. OR . leapf ) THEN 366 IF((.not.forward).OR. leapf ) THEN 367 ! Ehouarn: gather mass fluxes during backward Matsuno or LF step 361 368 CALL caladvtrac(q,pbaru,pbarv, 362 369 * p, masse, dq, teta, … … 497 504 ENDDO 498 505 ENDDO ! of DO l=1,llm 499 call friction(ucov,vcov,dtvr)500 506 501 if (planet_type.eq."giant") then 502 ! Intrinsic heat flux 503 ! Aymeric -- for giant planets 504 if (ihf .gt. 1.e-6) then 505 !print *, '**** INTRINSIC HEAT FLUX ****', ihf 506 DO ij=1,ip1jmp1 507 teta(ij,1) = teta(ij,1) 508 & + dtvr * aire(ij) * ihf / cpp / masse(ij,1) 509 ENDDO 510 !print *, '**** d teta ' 511 !print *, dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1) 512 endif 513 endif 507 if (planet_type.eq."giant") then 508 ! add an intrinsic heat flux at the base of the atmosphere 509 teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1) 510 endif 511 512 call friction(ucov,vcov,dtvr) 514 513 515 514 -
trunk/LMDZ.COMMON/libf/dyn3d/limit_netcdf.F90
r7 r270 1 1 ! 2 ! $Id: limit_netcdf.F90 1 441 2010-10-13 13:06:56Z emillour$2 ! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $ 3 3 !------------------------------------------------------------------------------- 4 4 ! … … 42 42 REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask 43 43 #ifndef CPP_EARTH 44 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'44 CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1) 45 45 #else 46 46 !------------------------------------------------------------------------------- … … 52 52 #include "indicesol.h" 53 53 54 !--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------55 LOGICAL, PARAMETER :: fracterre=.TRUE.56 57 54 !--- INPUT NETCDF FILES NAMES -------------------------------------------------- 58 55 CHARACTER(LEN=25) :: icefile, sstfile, dumstr 59 56 CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc ', & 60 57 famipsic='amipbc_sic_1x1.nc ', & 61 fclimsst='amipbc_sst_1x1_clim.nc ', &62 fclimsic='amipbc_sic_1x1_clim.nc ', &63 58 fcpldsst='cpl_atm_sst.nc ', & 64 59 fcpldsic='cpl_atm_sic.nc ', & 60 fhistsst='histmth_sst.nc ', & 61 fhistsic='histmth_sic.nc ', & 65 62 frugo ='Rugos.nc ', & 66 63 falbe ='Albedo.nc ' 67 64 CHARACTER(LEN=10) :: varname 68 65 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------ 69 66 REAL, DIMENSION(klon) :: fi_ice, verif … … 80 77 INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC 81 78 INTEGER :: NF90_FORMAT 82 LOGICAL :: lCPL !--- T: IPCC-IPSL cpl model output files83 79 INTEGER :: ndays !--- Depending on the output calendar 84 80 … … 104 100 105 101 !--- RUGOSITY TREATMENT -------------------------------------------------------- 106 WRITE(lunout,*) 'Traitement de la rugosite' 107 CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) 102 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite' 103 varname='RUGOS' 104 CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) 108 105 109 106 !--- OCEAN TREATMENT ----------------------------------------------------------- 110 PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.107 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique' 111 108 112 109 ! Input SIC file selection 113 icefile='fake' 114 IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic) 115 IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic) 116 IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic) 117 SELECT CASE(icefile) 118 CASE(famipsic); dumstr='Amip.' 119 CASE(fclimsic); dumstr='Amip climatologique.' 120 CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. 121 CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1) 122 END SELECT 110 ! Open file only to test if available 111 IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 112 icefile=TRIM(famipsic) 113 varname='sicbcs' 114 ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 115 icefile=TRIM(fcpldsic) 116 varname='SIICECOV' 117 ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 118 icefile=TRIM(fhistsic) 119 varname='pourc_sic' 120 ELSE 121 WRITE(lunout,*) 'ERROR! No sea-ice input file was found.' 122 WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic) 123 CALL abort_gcm('limit_netcdf','No sea-ice file was found',1) 124 END IF 123 125 ierr=NF90_CLOSE(nid) 124 WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)125 126 CALL get_2Dfield(icefile, 'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL)126 IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile) 127 128 CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice) 127 129 128 130 ALLOCATE(pctsrf_t(klon,nbsrf,ndays)) 129 131 DO k=1,ndays 130 fi_ice=phy_ice(:,k) 131 WHERE(fi_ice>=1.0 ) fi_ice=1.0 132 WHERE(fi_ice<EPSFRA) fi_ice=0.0 133 IF(fracterre) THEN 134 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil 135 pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice 136 IF(lCPL) THEN ! SIC=pICE*(1-LIC-TER) 137 pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter)) 138 ELSE ! SIC=pICE-LIC 132 fi_ice=phy_ice(:,k) 133 WHERE(fi_ice>=1.0 ) fi_ice=1.0 134 WHERE(fi_ice<EPSFRA) fi_ice=0.0 135 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil 136 pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice 137 IF (icefile==trim(fcpldsic)) THEN ! SIC=pICE*(1-LIC-TER) 138 pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter)) 139 ELSE IF (icefile==trim(fhistsic)) THEN ! SIC=pICE 140 pctsrf_t(:,is_sic,k)=fi_ice(:) 141 ELSE ! icefile==famipsic ! SIC=pICE-LIC 139 142 pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k) 140 141 142 143 END IF 144 WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0. 145 WHERE(1.0-zmasq<EPSFRA) 143 146 pctsrf_t(:,is_sic,k)=0.0 144 147 pctsrf_t(:,is_oce,k)=0.0 145 148 ELSEWHERE 146 149 WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq) 147 pctsrf_t(:,is_sic,k)=1.0-zmasq148 pctsrf_t(:,is_oce,k)=0.0150 pctsrf_t(:,is_sic,k)=1.0-zmasq 151 pctsrf_t(:,is_oce,k)=0.0 149 152 ELSEWHERE 150 pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)151 WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)152 pctsrf_t(:,is_oce,k)=0.0153 pctsrf_t(:,is_sic,k)=1.0-zmasq154 END WHERE153 pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k) 154 WHERE(pctsrf_t(:,is_oce,k)<EPSFRA) 155 pctsrf_t(:,is_oce,k)=0.0 156 pctsrf_t(:,is_sic,k)=1.0-zmasq 157 END WHERE 155 158 END WHERE 156 END WHERE 157 nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) 158 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 159 nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) 160 IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad 161 ELSE 162 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) 163 WHERE(NINT(pctsrf(:,is_ter))==1) 164 pctsrf_t(:,is_sic,k)=0. 165 pctsrf_t(:,is_oce,k)=0. 166 WHERE(fi_ice>=1.e-5) 167 pctsrf_t(:,is_lic,k)=fi_ice 168 pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k) 169 ELSEWHERE 170 pctsrf_t(:,is_lic,k)=0.0 171 END WHERE 172 ELSEWHERE 173 pctsrf_t(:,is_lic,k) = 0.0 174 WHERE(fi_ice>=1.e-5) 175 pctsrf_t(:,is_ter,k)=0.0 176 pctsrf_t(:,is_sic,k)=fi_ice 177 pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k) 178 ELSEWHERE 179 pctsrf_t(:,is_sic,k)=0.0 180 pctsrf_t(:,is_oce,k)=1.0 181 END WHERE 182 END WHERE 183 verif=sum(pctsrf_t(:,:,k),dim=2) 184 nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5) 185 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 186 END IF 159 END WHERE 160 nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) 161 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 162 nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) 163 IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad 187 164 END DO 188 165 DEALLOCATE(phy_ice) 189 166 190 167 !--- SST TREATMENT ------------------------------------------------------------- 191 WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.168 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst' 192 169 193 170 ! Input SST file selection 194 sstfile='fake' 195 IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst) 196 IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst) 197 IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst) 198 SELECT CASE(icefile) 199 CASE(famipsic); dumstr='Amip.' 200 CASE(fclimsic); dumstr='Amip climatologique.' 201 CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. 202 CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1) 203 END SELECT 171 ! Open file only to test if available 172 IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 173 sstfile=TRIM(famipsst) 174 varname='tosbcs' 175 ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 176 sstfile=TRIM(fcpldsst) 177 varname='SISUTESW' 178 ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 179 sstfile=TRIM(fhistsst) 180 varname='tsol_oce' 181 ELSE 182 WRITE(lunout,*) 'ERROR! No sst input file was found.' 183 WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst) 184 CALL abort_gcm('limit_netcdf','No sst file was found',1) 185 END IF 204 186 ierr=NF90_CLOSE(nid) 205 WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)206 207 CALL get_2Dfield( trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)187 IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile) 188 189 CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap) 208 190 209 191 !--- ALBEDO TREATMENT ---------------------------------------------------------- 210 WRITE(lunout,*) 'Traitement de l albedo' 211 CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb) 192 IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo' 193 varname='ALBEDO' 194 CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb) 212 195 213 196 !--- REFERENCE GROUND HEAT FLUX TREATMENT -------------------------------------- … … 215 198 216 199 !--- OUTPUT FILE WRITING ------------------------------------------------------- 217 WRITE(lunout,*) 'Ecriture du fichier limit'200 IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut' 218 201 219 202 !--- File creation … … 264 247 ierr=NF90_CLOSE(nid) 265 248 249 IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin' 250 266 251 DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug) 267 252 … … 276 261 !------------------------------------------------------------------------------- 277 262 ! 278 SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)263 SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask) 279 264 ! 280 265 !----------------------------------------------------------------------------- … … 304 289 ! Arguments: 305 290 CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name 291 CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name 306 292 CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB 307 293 LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels 308 294 INTEGER, INTENT(IN) :: ndays ! current year number of days 309 REAL, POINTER, DIMENSION(:, :) :: champo 295 REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t) 310 296 LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC) 311 297 REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask 312 LOGICAL, OPTIONAL, INTENT(IN) :: lCPL ! Coupled model flag (for ICE)313 298 !------------------------------------------------------------------------------ 314 299 ! Local variables: … … 316 301 INTEGER :: ncid, varid ! NetCDF identifiers 317 302 CHARACTER(LEN=30) :: dnam ! dimension name 318 CHARACTER(LEN=80) :: varname ! NetCDF variable name319 303 !--- dimensions 320 304 INTEGER, DIMENSION(4) :: dids ! NetCDF dimensions identifiers … … 331 315 !--- input files 332 316 CHARACTER(LEN=20) :: cal_in ! calendar 317 CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file 333 318 INTEGER :: ndays_in ! number of days 334 319 !--- misc … … 337 322 CHARACTER(LEN=25) :: title ! for messages 338 323 LOGICAL :: extrp ! flag for extrapolation 324 LOGICAL :: oldice ! flag for old way ice computation 339 325 REAL :: chmin, chmax 340 326 INTEGER ierr 341 327 integer n_extrap ! number of extrapolated points 342 328 logical skip 329 343 330 !------------------------------------------------------------------------------ 344 331 !---Variables depending on keyword 'mode' ------------------------------------- 345 332 NULLIFY(champo) 333 346 334 SELECT CASE(mode) 347 CASE('RUG'); varname='RUGOS';title='Rugosite'348 CASE('SIC'); varname='sicbcs'; title='Sea-ice'349 CASE('SST'); varname='tosbcs'; title='SST'350 CASE('ALB'); varname='ALBEDO'; title='Albedo'335 CASE('RUG'); title='Rugosite' 336 CASE('SIC'); title='Sea-ice' 337 CASE('SST'); title='SST' 338 CASE('ALB'); title='Albedo' 351 339 END SELECT 340 341 352 342 extrp=.FALSE. 343 oldice=.FALSE. 353 344 IF ( PRESENT(flag) ) THEN 354 345 IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 346 IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. 355 347 END IF 356 348 357 349 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ----------------------------- 350 IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam 358 351 ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid); CALL ncerr(ierr, fnam) 359 ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam)352 ierr=NF90_INQ_VARID(ncid, trim(varname), varid); CALL ncerr(ierr, fnam) 360 353 ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam) 354 355 !--- Read unit for sea-ice variable only 356 IF (mode=='SIC') THEN 357 ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic) 358 IF(ierr/=NF90_NOERR) THEN 359 IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value' 360 unit_sic='X' 361 ELSE 362 IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic 363 END IF 364 END IF 361 365 362 366 !--- Longitude … … 365 369 ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) 366 370 ierr=NF90_GET_VAR(ncid, varid, dlon_ini); CALL ncerr(ierr, fnam) 367 WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep371 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep 368 372 369 373 !--- Latitude … … 372 376 ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) 373 377 ierr=NF90_GET_VAR(ncid, varid, dlat_ini); CALL ncerr(ierr, fnam) 374 WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep378 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep 375 379 376 380 !--- Time (variable is not needed - it is rebuilt - but calendar is) … … 385 389 CASE('SIC', 'SST'); cal_in='gregorian' 386 390 END SELECT 387 WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &391 IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' & 388 392 // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.' 389 393 END IF 390 WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &394 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', & 391 395 cal_in 392 396 397 393 398 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION -------------------- 394 399 !--- Determining input file number of days, depending on calendar … … 398 403 !--- If input records are not monthly, time sampling has to be constant ! 399 404 timeyear=mid_months(anneeref, cal_in, lmdep) 400 IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' & 401 // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, & 402 ' enregistrements.' 405 IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), & 406 ' ne comportent pas 12, mais ', lmdep, ' enregistrements.' 403 407 404 408 !--- GETTING THE FIELD AND INTERPOLATING IT ---------------------------------- … … 406 410 IF(extrp) ALLOCATE(work(imdep, jmdep)) 407 411 408 WRITE(lunout, *) 409 WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, & 410 ' CHAMPS.' 412 IF (prt_level>5) WRITE(lunout, *) 413 IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.' 411 414 ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) 412 415 DO l=1, lmdep … … 419 422 work) 420 423 421 IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN 422 IF(l==1) THEN 423 WRITE(lunout, *) & 424 '-------------------------------------------------------------------------' 425 WRITE(lunout, *) & 426 'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$' 427 WRITE(lunout, *) & 428 '-------------------------------------------------------------------------' 424 IF(ibar .AND. .NOT.oldice) THEN 425 IF(l==1 .AND. prt_level>5) THEN 426 WRITE(lunout, *) '-------------------------------------------------------------------------' 427 WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$' 428 WRITE(lunout, *) '-------------------------------------------------------------------------' 429 429 END IF 430 430 IF(mode=='RUG') champ=LOG(champ) … … 453 453 454 454 !--- TIME INTERPOLATION ------------------------------------------------------ 455 WRITE(lunout, *) 456 WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' 457 WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear 458 WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays 455 IF (prt_level>5) THEN 456 WRITE(lunout, *) 457 WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' 458 WRITE(lunout, *)' Vecteur temps en entree: ', timeyear 459 WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays 460 END IF 461 459 462 ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays)) 460 463 skip = .false. … … 471 474 END DO 472 475 if (n_extrap /= 0) then 473 print *,"get_2Dfield pchfe_95: n_extrap = ", n_extrap476 WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap 474 477 end if 475 478 champan(iip1, :, :)=champan(1, :, :) … … 479 482 DO j=1, jjp1 480 483 CALL minmax(iip1, champan(1, j, 10), chmin, chmax) 481 WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j484 IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j 482 485 END DO 483 486 484 487 !--- SPECIAL FILTER FOR SST: SST>271.38 -------------------------------------- 485 488 IF(mode=='SST') THEN 486 WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'489 IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38' 487 490 WHERE(champan<271.38) champan=271.38 488 491 END IF … … 490 493 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ------------------------------------- 491 494 IF(mode=='SIC') THEN 492 WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0' 493 IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100. 495 IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0' 496 497 IF (unit_sic=='1') THEN 498 ! Nothing to be done for sea-ice field is already in fraction of 1 499 ! This is the case for sea-ice in file cpl_atm_sic.nc 500 IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1' 501 ELSE 502 ! Convert sea ice from percentage to fraction of 1 503 IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' 504 champan(:, :, :)=champan(:, :, :)/100. 505 END IF 506 494 507 champan(iip1, :, :)=champan(1, :, :) 495 508 WHERE(champan>1.0) champan=1.0 496 509 WHERE(champan<0.0) champan=0.0 497 510 END IF 498 511 499 512 !--- DYNAMICAL TO PHYSICAL GRID ---------------------------------------------- -
trunk/LMDZ.COMMON/libf/dyn3d/vlsplt.F
r109 r270 1 ! 2 ! $Header$ 3 ! 4 c 1 c 2 c $Id: vlsplt.F 1550 2011-07-05 09:44:55Z lguez $ 5 3 c 6 4 … … 131 129 IMPLICIT NONE 132 130 c 133 #include "dimensions.h" 134 #include "paramet.h" 135 #include "logic.h" 136 #include "comvert.h" 137 #include "comconst.h" 131 include "dimensions.h" 132 include "paramet.h" 133 include "logic.h" 134 include "comvert.h" 135 include "comconst.h" 136 include "iniprint.h" 138 137 c 139 138 c … … 353 352 354 353 IF(n0.gt.0) THEN 355 PRINT*,'Nombre de points pour lesquels on advect plus que le' 354 if (prt_level > 2) PRINT *, 355 $ 'Nombre de points pour lesquels on advect plus que le' 356 356 & ,'contenu de la maille : ',n0 357 357 -
trunk/LMDZ.COMMON/libf/dyn3d/write_paramLMDZ_dyn.h
r1 r270 51 51 . zx_tmp_2d,iip1*jjp1,ndex2d) 52 52 c 53 zx_tmp_2d(1:iip1,1:jjp1)=REAL( idissip)54 CALL histwrite(nid_ctesGCM, " idissip", itau_w,53 zx_tmp_2d(1:iip1,1:jjp1)=REAL(dissip_period) 54 CALL histwrite(nid_ctesGCM, "dissip_period", itau_w, 55 55 . zx_tmp_2d,iip1*jjp1,ndex2d) 56 56 c
Note: See TracChangeset
for help on using the changeset viewer.
