Changeset 2720 for LMDZ5/branches/testing/libf/phylmd/dyn1d
- Timestamp:
- Nov 30, 2016, 1:28:41 PM (8 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 7 edited
- 1 copied
-
. (modified) (1 prop)
-
libf/phylmd/dyn1d/1DUTILS.h (modified) (21 diffs)
-
libf/phylmd/dyn1d/1D_decl_cases.h (modified) (7 diffs)
-
libf/phylmd/dyn1d/1D_interp_cases.h (modified) (3 diffs)
-
libf/phylmd/dyn1d/1D_read_forc_cases.h (modified) (6 diffs)
-
libf/phylmd/dyn1d/compar1d.h (modified) (4 diffs)
-
libf/phylmd/dyn1d/lmdz1d.F90 (modified) (18 diffs)
-
libf/phylmd/dyn1d/mod_1D_cases_read2.F90 (copied) (copied from LMDZ5/trunk/libf/phylmd/dyn1d/mod_1D_cases_read2.F90)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 2665-2668,2670-2674,2677-2681,2683-2684,2686,2690-2719
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/phylmd/dyn1d/1DUTILS.h
r2641 r2720 55 55 56 56 !Config Key = prt_level 57 !Config Desc = niveau d'impressions de d ébogage57 !Config Desc = niveau d'impressions de debogage 58 58 !Config Def = 0 59 !Config Help = Niveau d'impression pour le d ébogage59 !Config Help = Niveau d'impression pour le debogage 60 60 !Config (0 = minimum d'impression) 61 61 ! prt_level = 0 … … 118 118 ! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s 119 119 ! Radiation to be switched off 120 ! > 100 ==> forcing_case = .true. or forcing_case2 = .true. 121 ! initial profiles from case.nc file 120 122 ! 121 123 forcing_type = 0 … … 134 136 ENDIF 135 137 136 !Param ètres de forçage138 !Parametres de forcage 137 139 !Config Key = tend_t 138 140 !Config Desc = forcage ou non par advection de T … … 303 305 CALL getin('rugos',rugos) 304 306 307 !Config Key = rugosh 308 !Config Desc = coefficient de frottement 309 !Config Def = rugos 310 !Config Help = calcul du Cdrag 311 rugosh = rugos 312 CALL getin('rugosh',rugosh) 313 314 315 316 !Config Key = snowmass 317 !Config Desc = mass de neige de la surface en kg/m2 318 !Config Def = 0.0000 319 !Config Help = snowmass 320 snowmass = 0.0000 321 CALL getin('snowmass',snowmass) 322 305 323 !Config Key = wtsurf et wqsurf 306 324 !Config Desc = ??? … … 342 360 !Config Key = zpicinp 343 361 !Config Desc = denivellation orographie 344 !Config Def = 300.362 !Config Def = 0. 345 363 !Config Help = input brise 346 zpicinp = 300.364 zpicinp = 0. 347 365 CALL getin('zpicinp',zpicinp) 348 366 !Config key = nudge_tsoil … … 378 396 CALL getin('tau_soil_nudge',tau_soil_nudge) 379 397 398 !---------------------------------------------------------- 399 ! Param??tres de for??age pour les forcages communs: 400 ! Pour les forcages communs: ces entiers valent 0 ou 1 401 ! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale 402 ! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale 403 ! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv) 404 ! forcages en omega, w, vent geostrophique ou ustar 405 ! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging 406 !---------------------------------------------------------- 407 408 !Config Key = tadv 409 !Config Desc = forcage ou non par advection totale de T 410 !Config Def = false 411 !Config Help = forcage ou non par advection totale de T 412 tadv =0 413 CALL getin('tadv',tadv) 414 415 !Config Key = tadvv 416 !Config Desc = forcage ou non par advection verticale de T 417 !Config Def = false 418 !Config Help = forcage ou non par advection verticale de T 419 tadvv =0 420 CALL getin('tadvv',tadvv) 421 422 !Config Key = tadvh 423 !Config Desc = forcage ou non par advection horizontale de T 424 !Config Def = false 425 !Config Help = forcage ou non par advection horizontale de T 426 tadvh =0 427 CALL getin('tadvh',tadvh) 428 429 !Config Key = thadv 430 !Config Desc = forcage ou non par advection totale de Theta 431 !Config Def = false 432 !Config Help = forcage ou non par advection totale de Theta 433 thadv =0 434 CALL getin('thadv',thadv) 435 436 !Config Key = thadvv 437 !Config Desc = forcage ou non par advection verticale de Theta 438 !Config Def = false 439 !Config Help = forcage ou non par advection verticale de Theta 440 thadvv =0 441 CALL getin('thadvv',thadvv) 442 443 !Config Key = thadvh 444 !Config Desc = forcage ou non par advection horizontale de Theta 445 !Config Def = false 446 !Config Help = forcage ou non par advection horizontale de Theta 447 thadvh =0 448 CALL getin('thadvh',thadvh) 449 450 !Config Key = qadv 451 !Config Desc = forcage ou non par advection totale de Q 452 !Config Def = false 453 !Config Help = forcage ou non par advection totale de Q 454 qadv =0 455 CALL getin('qadv',qadv) 456 457 !Config Key = qadvv 458 !Config Desc = forcage ou non par advection verticale de Q 459 !Config Def = false 460 !Config Help = forcage ou non par advection verticale de Q 461 qadvv =0 462 CALL getin('qadvv',qadvv) 463 464 !Config Key = qadvh 465 !Config Desc = forcage ou non par advection horizontale de Q 466 !Config Def = false 467 !Config Help = forcage ou non par advection horizontale de Q 468 qadvh =0 469 CALL getin('qadvh',qadvh) 470 471 !Config Key = trad 472 !Config Desc = forcage ou non par tendance radiative 473 !Config Def = false 474 !Config Help = forcage ou non par tendance radiative 475 trad =0 476 CALL getin('trad',trad) 477 478 !Config Key = forc_omega 479 !Config Desc = forcage ou non par omega 480 !Config Def = false 481 !Config Help = forcage ou non par omega 482 forc_omega =0 483 CALL getin('forc_omega',forc_omega) 484 485 !Config Key = forc_w 486 !Config Desc = forcage ou non par w 487 !Config Def = false 488 !Config Help = forcage ou non par w 489 forc_w =0 490 CALL getin('forc_w',forc_w) 491 492 !Config Key = forc_geo 493 !Config Desc = forcage ou non par geo 494 !Config Def = false 495 !Config Help = forcage ou non par geo 496 forc_geo =0 497 CALL getin('forc_geo',forc_geo) 498 499 ! Meme chose que ok_precr_ust 500 !Config Key = forc_ustar 501 !Config Desc = forcage ou non par ustar 502 !Config Def = false 503 !Config Help = forcage ou non par ustar 504 forc_ustar =0 505 CALL getin('forc_ustar',forc_ustar) 506 IF (forc_ustar .EQ. 1) ok_prescr_ust=.true. 507 508 !Config Key = nudging_u 509 !Config Desc = forcage ou non par nudging sur u 510 !Config Def = false 511 !Config Help = forcage ou non par nudging sur u 512 nudging_u =0 513 CALL getin('nudging_u',nudging_u) 514 515 !Config Key = nudging_v 516 !Config Desc = forcage ou non par nudging sur v 517 !Config Def = false 518 !Config Help = forcage ou non par nudging sur v 519 nudging_v =0 520 CALL getin('nudging_v',nudging_v) 521 522 !Config Key = nudging_w 523 !Config Desc = forcage ou non par nudging sur w 524 !Config Def = false 525 !Config Help = forcage ou non par nudging sur w 526 nudging_w =0 527 CALL getin('nudging_w',nudging_w) 528 529 !Config Key = nudging_q 530 !Config Desc = forcage ou non par nudging sur q 531 !Config Def = false 532 !Config Help = forcage ou non par nudging sur q 533 nudging_q =0 534 CALL getin('nudging_q',nudging_q) 535 536 !Config Key = nudging_t 537 !Config Desc = forcage ou non par nudging sur t 538 !Config Def = false 539 !Config Help = forcage ou non par nudging sur t 540 nudging_t =0 541 CALL getin('nudging_t',nudging_t) 380 542 381 543 … … 395 557 write(lunout,*)' zsurf = ', zsurf 396 558 write(lunout,*)' rugos = ', rugos 559 write(lunout,*)' snowmass=', snowmass 397 560 write(lunout,*)' wtsurf = ', wtsurf 398 561 write(lunout,*)' wqsurf = ', wqsurf … … 406 569 write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge 407 570 write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge 571 write(lunout,*)' tadv = ', tadv 572 write(lunout,*)' tadvv = ', tadvv 573 write(lunout,*)' tadvh = ', tadvh 574 write(lunout,*)' thadv = ', thadv 575 write(lunout,*)' thadvv = ', thadvv 576 write(lunout,*)' thadvh = ', thadvh 577 write(lunout,*)' qadv = ', qadv 578 write(lunout,*)' qadvv = ', qadvv 579 write(lunout,*)' qadvh = ', qadvh 580 write(lunout,*)' trad = ', trad 581 write(lunout,*)' forc_omega = ', forc_omega 582 write(lunout,*)' forc_w = ', forc_w 583 write(lunout,*)' forc_geo = ', forc_geo 584 write(lunout,*)' forc_ustar = ', forc_ustar 585 write(lunout,*)' nudging_u = ', nudging_u 586 write(lunout,*)' nudging_v = ', nudging_v 587 write(lunout,*)' nudging_t = ', nudging_t 588 write(lunout,*)' nudging_q = ', nudging_q 408 589 IF (forcing_type .eq.40) THEN 409 590 write(lunout,*) '--- Forcing type GCSS Old --- with:' … … 1106 1287 !---------------------------------------------------------------------- 1107 1288 ! Calcul de l'advection verticale (ascendance et subsidence) de 1108 ! temp érature et d'humidité. Hypothèse : ce qui rentre de l'extérieur1109 ! a les m êmes caractéristiques que l'air de la colonne 1D (WTG) ou1289 ! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur 1290 ! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 1110 1291 ! sans WTG rajouter une advection horizontale 1111 1292 !---------------------------------------------------------------------- … … 1180 1361 !---------------------------------------------------------------------- 1181 1362 ! Calcul de l'advection verticale (ascendance et subsidence) de 1182 ! temp érature et d'humidité. Hypothèse : ce qui rentre de l'extérieur1183 ! a les m êmes caractéristiques que l'air de la colonne 1D (WTG) ou1363 ! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur 1364 ! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 1184 1365 ! sans WTG rajouter une advection horizontale 1185 1366 !---------------------------------------------------------------------- … … 2934 3115 endif 2935 3116 if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then 2936 print*,'TOGA-COARE a d ébutéle 1er Nov 1992 (jour julien=306)'3117 print*,'TOGA-COARE a debute le 1er Nov 1992 (jour julien=306)' 2937 3118 print*,'Changer dayref dans run.def' 2938 3119 stop … … 3141 3322 3142 3323 !====================================================================== 3324 SUBROUTINE interp_gabls4_time(day,day1,annee_ref & 3325 & ,year_ini_gabls4,day_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4 & 3326 & ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4 & 3327 & ,ug_prof,vg_prof,ht_prof,hq_prof,tg_prof) 3328 implicit none 3329 3330 !--------------------------------------------------------------------------------------- 3331 ! Time interpolation of a 2D field to the timestep corresponding to day 3332 ! 3333 ! day: current julian day 3334 ! day1: first day of the simulation 3335 ! nt_gabls4: total nb of data in the forcing (e.g. 37 for gabls4) 3336 ! dt_gabls4: total time interval (in sec) between 2 forcing data (e.g. 60min. for gabls4) 3337 !--------------------------------------------------------------------------------------- 3338 3339 #include "compar1d.h" 3340 3341 ! inputs: 3342 integer annee_ref 3343 integer nt_gabls4,nlev_gabls4 3344 integer year_ini_gabls4 3345 real day, day1,day_ini_gabls4,dt_gabls4 3346 real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4) 3347 real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4) 3348 real tg_gabls4(nt_gabls4), tg_prof 3349 ! outputs: 3350 real ug_prof(nlev_gabls4),vg_prof(nlev_gabls4) 3351 real ht_prof(nlev_gabls4),hq_prof(nlev_gabls4) 3352 ! local: 3353 integer it_gabls41, it_gabls42,k 3354 real timeit,time_gabls41,time_gabls42,frac 3355 3356 3357 3358 ! Check that initial day of the simulation consistent with gabls4 period: 3359 if (forcing_type.eq.8 ) then 3360 print *,'annee_ref=',annee_ref 3361 print *,'day1=',day1 3362 print *,'day_ini_gabls4=',day_ini_gabls4 3363 if (annee_ref.ne.2009) then 3364 print*,'Pour gabls4, annee_ref doit etre 2009' 3365 print*,'Changer annee_ref dans run.def' 3366 stop 3367 endif 3368 if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4) then 3369 print*,'gabls4 a debute le 11 dec 2009 (jour julien=345)' 3370 print*,'Changer dayref dans run.def',day1,day_ini_gabls4 3371 stop 3372 endif 3373 if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4+2) then 3374 print*,'gabls4 a fini le 12 dec 2009 (jour julien=346)' 3375 print*,'Changer dayref ou nday dans run.def',day1,day_ini_gabls4 3376 stop 3377 endif 3378 endif 3379 3380 timeit=(day-day_ini_gabls4)*86400 3381 print *,'day,day_ini_gabls4=',day,day_ini_gabls4 3382 print *,'nt_gabls4,dt,timeit=',nt_gabls4,dt_gabls4,timeit 3383 3384 ! Determine the closest observation times: 3385 it_gabls41=INT(timeit/dt_gabls4)+1 3386 it_gabls42=it_gabls41 + 1 3387 time_gabls41=(it_gabls41-1)*dt_gabls4 3388 time_gabls42=(it_gabls42-1)*dt_gabls4 3389 3390 if (it_gabls41 .ge. nt_gabls4) then 3391 write(*,*) 'PB-stop: day, it_gabls41, it_gabls42, timeit: ',day,it_gabls41,it_gabls42,timeit/86400. 3392 stop 3393 endif 3394 3395 ! time interpolation: 3396 frac=(time_gabls42-timeit)/(time_gabls42-time_gabls41) 3397 frac=max(frac,0.0) 3398 3399 3400 do k=1,nlev_gabls4 3401 ug_prof(k) = ug_gabls4(k,it_gabls42)-frac*(ug_gabls4(k,it_gabls42)-ug_gabls4(k,it_gabls41)) 3402 vg_prof(k) = vg_gabls4(k,it_gabls42)-frac*(vg_gabls4(k,it_gabls42)-vg_gabls4(k,it_gabls41)) 3403 ht_prof(k) = ht_gabls4(k,it_gabls42)-frac*(ht_gabls4(k,it_gabls42)-ht_gabls4(k,it_gabls41)) 3404 hq_prof(k) = hq_gabls4(k,it_gabls42)-frac*(hq_gabls4(k,it_gabls42)-hq_gabls4(k,it_gabls41)) 3405 enddo 3406 tg_prof=tg_gabls4(it_gabls42)-frac*(tg_gabls4(it_gabls42)-tg_gabls4(it_gabls41)) 3407 return 3408 END 3409 3410 !====================================================================== 3143 3411 SUBROUTINE interp_armcu_time(day,day1,annee_ref & 3144 3412 & ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu & … … 3679 3947 !===================================================================== 3680 3948 subroutine read_dice(fich_dice,nlevel,ntime & 3681 & ,zz,pres,t h,qv,u,v,o3 &3949 & ,zz,pres,t,qv,u,v,o3 & 3682 3950 & ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg & 3683 3951 & ,hadvt,hadvq,hadvu,hadvv,w,omega) … … 3689 3957 3690 3958 #include "netcdf.inc" 3959 #include "YOMCST.h" 3691 3960 3692 3961 integer ntime,nlevel … … 3696 3965 real*8 zz(nlevel) 3697 3966 3698 real*8 th(nlevel),pres(nlevel) 3967 real*8 th(nlevel),pres(nlevel),t(nlevel) 3699 3968 real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel) 3700 3969 real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime) … … 3702 3971 real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime) 3703 3972 real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime) 3973 real*8 pzero 3704 3974 3705 3975 integer nid, ierr … … 3708 3978 integer var3didin(nbvar3d) 3709 3979 3980 pzero=100000. 3710 3981 ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid) 3711 3982 if (ierr.NE.NF_NOERR) then … … 3882 4153 endif 3883 4154 ! write(*,*)'lecture th ok',th 4155 do k=1,nlevel 4156 t(k)=th(k)*(pres(k)/pzero)**rkappa 4157 enddo 3884 4158 3885 4159 #ifdef NC_DOUBLE … … 4095 4369 end subroutine read_dice 4096 4370 !===================================================================== 4371 subroutine read_gabls4(fich_gabls4,nlevel,ntime,nsol & 4372 & ,zz,depth_sn,ug,vg,pf,th,t,qv,u,v,hadvt,hadvq,tg,tsnow,snow_dens) 4373 4374 !program reading initial profils and forcings of the Gabls4 case study 4375 4376 4377 implicit none 4378 4379 #include "netcdf.inc" 4380 4381 integer ntime,nlevel,nsol 4382 integer l,k 4383 character*80 :: fich_gabls4 4384 real*8 time(ntime) 4385 4386 ! ATTENTION: visiblement quand on lit gabls4_driver.nc on recupere les donnees 4387 ! dans un ordre inverse par rapport a la convention LMDZ 4388 ! ==> il faut tout inverser (MPL 20141024) 4389 ! les variables indexees "_i" sont celles qui sont lues dans gabls4_driver.nc 4390 real*8 zz_i(nlevel),th_i(nlevel),pf_i(nlevel),t_i(nlevel) 4391 real*8 qv_i(nlevel),u_i(nlevel),v_i(nlevel),ug_i(nlevel,ntime),vg_i(nlevel,ntime) 4392 real*8 hadvt_i(nlevel,ntime),hadvq_i(nlevel,ntime) 4393 4394 real*8 zz(nlevel),th(nlevel),pf(nlevel),t(nlevel) 4395 real*8 qv(nlevel),u(nlevel),v(nlevel),ug(nlevel,ntime),vg(nlevel,ntime) 4396 real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime) 4397 4398 real*8 depth_sn(nsol),tsnow(nsol),snow_dens(nsol) 4399 real*8 tg(ntime) 4400 integer nid, ierr 4401 integer nbvar3d 4402 parameter(nbvar3d=30) 4403 integer var3didin(nbvar3d) 4404 4405 ierr = NF_OPEN(fich_gabls4,NF_NOWRITE,nid) 4406 if (ierr.NE.NF_NOERR) then 4407 write(*,*) 'ERROR: Pb opening forcings nc file ' 4408 write(*,*) NF_STRERROR(ierr) 4409 stop "" 4410 endif 4411 4412 4413 ierr=NF_INQ_VARID(nid,"height",var3didin(1)) 4414 if(ierr/=NF_NOERR) then 4415 write(*,*) NF_STRERROR(ierr) 4416 stop 'height' 4417 endif 4418 4419 ierr=NF_INQ_VARID(nid,"depth_sn",var3didin(2)) 4420 if(ierr/=NF_NOERR) then 4421 write(*,*) NF_STRERROR(ierr) 4422 stop 'depth_sn' 4423 endif 4424 4425 ierr=NF_INQ_VARID(nid,"Ug",var3didin(3)) 4426 if(ierr/=NF_NOERR) then 4427 write(*,*) NF_STRERROR(ierr) 4428 stop 'Ug' 4429 endif 4430 4431 ierr=NF_INQ_VARID(nid,"Vg",var3didin(4)) 4432 if(ierr/=NF_NOERR) then 4433 write(*,*) NF_STRERROR(ierr) 4434 stop 'Vg' 4435 endif 4436 ierr=NF_INQ_VARID(nid,"pf",var3didin(5)) 4437 if(ierr/=NF_NOERR) then 4438 write(*,*) NF_STRERROR(ierr) 4439 stop 'pf' 4440 endif 4441 4442 ierr=NF_INQ_VARID(nid,"theta",var3didin(6)) 4443 if(ierr/=NF_NOERR) then 4444 write(*,*) NF_STRERROR(ierr) 4445 stop 'theta' 4446 endif 4447 4448 ierr=NF_INQ_VARID(nid,"tempe",var3didin(7)) 4449 if(ierr/=NF_NOERR) then 4450 write(*,*) NF_STRERROR(ierr) 4451 stop 'tempe' 4452 endif 4453 4454 ierr=NF_INQ_VARID(nid,"qv",var3didin(8)) 4455 if(ierr/=NF_NOERR) then 4456 write(*,*) NF_STRERROR(ierr) 4457 stop 'qv' 4458 endif 4459 4460 ierr=NF_INQ_VARID(nid,"u",var3didin(9)) 4461 if(ierr/=NF_NOERR) then 4462 write(*,*) NF_STRERROR(ierr) 4463 stop 'u' 4464 endif 4465 4466 ierr=NF_INQ_VARID(nid,"v",var3didin(10)) 4467 if(ierr/=NF_NOERR) then 4468 write(*,*) NF_STRERROR(ierr) 4469 stop 'v' 4470 endif 4471 4472 ierr=NF_INQ_VARID(nid,"hadvT",var3didin(11)) 4473 if(ierr/=NF_NOERR) then 4474 write(*,*) NF_STRERROR(ierr) 4475 stop 'hadvt' 4476 endif 4477 4478 ierr=NF_INQ_VARID(nid,"hadvQ",var3didin(12)) 4479 if(ierr/=NF_NOERR) then 4480 write(*,*) NF_STRERROR(ierr) 4481 stop 'hadvq' 4482 endif 4483 4484 ierr=NF_INQ_VARID(nid,"Tsnow",var3didin(14)) 4485 if(ierr/=NF_NOERR) then 4486 write(*,*) NF_STRERROR(ierr) 4487 stop 'tsnow' 4488 endif 4489 4490 ierr=NF_INQ_VARID(nid,"snow_density",var3didin(15)) 4491 if(ierr/=NF_NOERR) then 4492 write(*,*) NF_STRERROR(ierr) 4493 stop 'snow_density' 4494 endif 4495 4496 ierr=NF_INQ_VARID(nid,"Tg",var3didin(16)) 4497 if(ierr/=NF_NOERR) then 4498 write(*,*) NF_STRERROR(ierr) 4499 stop 'Tg' 4500 endif 4501 4502 4503 !dimensions lecture 4504 ! call catchaxis(nid,ntime,nlevel,time,z,ierr) 4505 4506 #ifdef NC_DOUBLE 4507 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz_i) 4508 #else 4509 ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz_i) 4510 #endif 4511 if(ierr/=NF_NOERR) then 4512 write(*,*) NF_STRERROR(ierr) 4513 stop "getvarup" 4514 endif 4515 4516 #ifdef NC_DOUBLE 4517 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),depth_sn) 4518 #else 4519 ierr = NF_GET_VAR_REAL(nid,var3didin(2),depth_sn) 4520 #endif 4521 if(ierr/=NF_NOERR) then 4522 write(*,*) NF_STRERROR(ierr) 4523 stop "getvarup" 4524 endif 4525 4526 #ifdef NC_DOUBLE 4527 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),ug_i) 4528 #else 4529 ierr = NF_GET_VAR_REAL(nid,var3didin(3),ug_i) 4530 #endif 4531 if(ierr/=NF_NOERR) then 4532 write(*,*) NF_STRERROR(ierr) 4533 stop "getvarup" 4534 endif 4535 4536 #ifdef NC_DOUBLE 4537 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),vg_i) 4538 #else 4539 ierr = NF_GET_VAR_REAL(nid,var3didin(4),vg_i) 4540 #endif 4541 if(ierr/=NF_NOERR) then 4542 write(*,*) NF_STRERROR(ierr) 4543 stop "getvarup" 4544 endif 4545 4546 #ifdef NC_DOUBLE 4547 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),pf_i) 4548 #else 4549 ierr = NF_GET_VAR_REAL(nid,var3didin(5),pf_i) 4550 #endif 4551 if(ierr/=NF_NOERR) then 4552 write(*,*) NF_STRERROR(ierr) 4553 stop "getvarup" 4554 endif 4555 4556 #ifdef NC_DOUBLE 4557 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),th_i) 4558 #else 4559 ierr = NF_GET_VAR_REAL(nid,var3didin(6),th_i) 4560 #endif 4561 if(ierr/=NF_NOERR) then 4562 write(*,*) NF_STRERROR(ierr) 4563 stop "getvarup" 4564 endif 4565 4566 #ifdef NC_DOUBLE 4567 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),t_i) 4568 #else 4569 ierr = NF_GET_VAR_REAL(nid,var3didin(7),t_i) 4570 #endif 4571 if(ierr/=NF_NOERR) then 4572 write(*,*) NF_STRERROR(ierr) 4573 stop "getvarup" 4574 endif 4575 4576 #ifdef NC_DOUBLE 4577 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),qv_i) 4578 #else 4579 ierr = NF_GET_VAR_REAL(nid,var3didin(8),qv_i) 4580 #endif 4581 if(ierr/=NF_NOERR) then 4582 write(*,*) NF_STRERROR(ierr) 4583 stop "getvarup" 4584 endif 4585 4586 #ifdef NC_DOUBLE 4587 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),u_i) 4588 #else 4589 ierr = NF_GET_VAR_REAL(nid,var3didin(9),u_i) 4590 #endif 4591 if(ierr/=NF_NOERR) then 4592 write(*,*) NF_STRERROR(ierr) 4593 stop "getvarup" 4594 endif 4595 4596 #ifdef NC_DOUBLE 4597 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),v_i) 4598 #else 4599 ierr = NF_GET_VAR_REAL(nid,var3didin(10),v_i) 4600 #endif 4601 if(ierr/=NF_NOERR) then 4602 write(*,*) NF_STRERROR(ierr) 4603 stop "getvarup" 4604 endif 4605 4606 #ifdef NC_DOUBLE 4607 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),hadvt_i) 4608 #else 4609 ierr = NF_GET_VAR_REAL(nid,var3didin(11),hadvt_i) 4610 #endif 4611 if(ierr/=NF_NOERR) then 4612 write(*,*) NF_STRERROR(ierr) 4613 stop "getvarup" 4614 endif 4615 4616 #ifdef NC_DOUBLE 4617 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),hadvq_i) 4618 #else 4619 ierr = NF_GET_VAR_REAL(nid,var3didin(12),hadvq_i) 4620 #endif 4621 if(ierr/=NF_NOERR) then 4622 write(*,*) NF_STRERROR(ierr) 4623 stop "getvarup" 4624 endif 4625 4626 #ifdef NC_DOUBLE 4627 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),tsnow) 4628 #else 4629 ierr = NF_GET_VAR_REAL(nid,var3didin(14),tsnow) 4630 #endif 4631 if(ierr/=NF_NOERR) then 4632 write(*,*) NF_STRERROR(ierr) 4633 stop "getvarup" 4634 endif 4635 4636 #ifdef NC_DOUBLE 4637 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),snow_dens) 4638 #else 4639 ierr = NF_GET_VAR_REAL(nid,var3didin(15),snow_dens) 4640 #endif 4641 if(ierr/=NF_NOERR) then 4642 write(*,*) NF_STRERROR(ierr) 4643 stop "getvarup" 4644 endif 4645 4646 #ifdef NC_DOUBLE 4647 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),tg) 4648 #else 4649 ierr = NF_GET_VAR_REAL(nid,var3didin(16),tg) 4650 #endif 4651 if(ierr/=NF_NOERR) then 4652 write(*,*) NF_STRERROR(ierr) 4653 stop "getvarup" 4654 endif 4655 4656 ! On remet les variables lues dans le bon ordre des niveaux (MPL 20141024) 4657 do k=1,nlevel 4658 zz(k)=zz_i(nlevel+1-k) 4659 ug(k,:)=ug_i(nlevel+1-k,:) 4660 vg(k,:)=vg_i(nlevel+1-k,:) 4661 pf(k)=pf_i(nlevel+1-k) 4662 print *,'pf=',pf(k) 4663 th(k)=th_i(nlevel+1-k) 4664 t(k)=t_i(nlevel+1-k) 4665 qv(k)=qv_i(nlevel+1-k) 4666 u(k)=u_i(nlevel+1-k) 4667 v(k)=v_i(nlevel+1-k) 4668 hadvt(k,:)=hadvt_i(nlevel+1-k,:) 4669 hadvq(k,:)=hadvq_i(nlevel+1-k,:) 4670 enddo 4671 return 4672 end subroutine read_gabls4 4673 !===================================================================== 4674 4097 4675 ! Reads CIRC input files 4098 4676 … … 4390 4968 ! 4391 4969 ! Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new) 4392 ! qui n' était pas correcte.4970 ! qui n'etait pas correcte. 4393 4971 ! 4394 4972 IF (tnew.LT.RTT) THEN … … 4465 5043 END 4466 5044 5045 !===================================================================== 5046 SUBROUTINE interp2_case_vertical(play,nlev_cas,plev_prof_cas & 5047 & ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas & 5048 & ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas & 5049 & ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas & 5050 & ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas & 5051 & ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & 5052 & ,dth_prof_cas,hth_prof_cas,vth_prof_cas & 5053 ! 5054 & ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas & 5055 & ,qv_mod_cas,ql_mod_cas,qi_mod_cas,u_mod_cas,v_mod_cas & 5056 & ,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas & 5057 & ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas & 5058 & ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas & 5059 & ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc) 5060 5061 implicit none 5062 5063 #include "dimensions.h" 5064 5065 !------------------------------------------------------------------------- 5066 ! Vertical interpolation of generic case forcing data onto mod_casel levels 5067 !------------------------------------------------------------------------- 5068 5069 integer nlevmax 5070 parameter (nlevmax=41) 5071 integer nlev_cas,mxcalc 5072 ! real play(llm), plev_prof(nlevmax) 5073 ! real t_prof(nlevmax),q_prof(nlevmax) 5074 ! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax) 5075 ! real ht_prof(nlevmax),vt_prof(nlevmax) 5076 ! real hq_prof(nlevmax),vq_prof(nlevmax) 5077 5078 real play(llm), plev_prof_cas(nlev_cas) 5079 real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas) 5080 real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas) 5081 real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas) 5082 real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas) 5083 real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas) 5084 real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas) 5085 real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas) 5086 real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas) 5087 real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas) 5088 5089 real t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm) 5090 real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm) 5091 real u_mod_cas(llm),v_mod_cas(llm) 5092 real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm) 5093 real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm) 5094 real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm) 5095 real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm) 5096 real dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm) 5097 real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm) 5098 5099 integer l,k,k1,k2 5100 real frac,frac1,frac2,fact 5101 5102 do l = 1, llm 5103 print *,'debut interp2, play=',l,play(l) 5104 enddo 5105 ! do l = 1, nlev_cas 5106 ! print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l) 5107 ! enddo 5108 5109 do l = 1, llm 5110 5111 if (play(l).ge.plev_prof_cas(nlev_cas)) then 5112 5113 mxcalc=l 5114 print *,'debut interp2, mxcalc=',mxcalc 5115 k1=0 5116 k2=0 5117 5118 if (play(l).le.plev_prof_cas(1)) then 5119 5120 do k = 1, nlev_cas-1 5121 if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then 5122 k1=k 5123 k2=k+1 5124 endif 5125 enddo 5126 5127 if (k1.eq.0 .or. k2.eq.0) then 5128 write(*,*) 'PB! k1, k2 = ',k1,k2 5129 write(*,*) 'l,play(l) = ',l,play(l)/100 5130 do k = 1, nlev_cas-1 5131 write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100 5132 enddo 5133 endif 5134 5135 frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1)) 5136 t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1)) 5137 theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1)) 5138 thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1)) 5139 thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1)) 5140 qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1)) 5141 ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1)) 5142 qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1)) 5143 u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1)) 5144 v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1)) 5145 ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1)) 5146 vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1)) 5147 w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1)) 5148 omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1)) 5149 du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1)) 5150 hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1)) 5151 vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1)) 5152 dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1)) 5153 hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1)) 5154 vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1)) 5155 dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1)) 5156 ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1)) 5157 vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1)) 5158 dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1)) 5159 hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1)) 5160 vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1)) 5161 dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1)) 5162 hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1)) 5163 vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1)) 5164 5165 else !play>plev_prof_cas(1) 5166 5167 k1=1 5168 k2=2 5169 print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2) 5170 frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2)) 5171 frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2)) 5172 t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2) 5173 theta_mod_cas(l)= frac1*th_prof_cas(k1) - frac2*th_prof_cas(k2) 5174 thv_mod_cas(l)= frac1*thv_prof_cas(k1) - frac2*thv_prof_cas(k2) 5175 thl_mod_cas(l)= frac1*thl_prof_cas(k1) - frac2*thl_prof_cas(k2) 5176 qv_mod_cas(l)= frac1*qv_prof_cas(k1) - frac2*qv_prof_cas(k2) 5177 ql_mod_cas(l)= frac1*ql_prof_cas(k1) - frac2*ql_prof_cas(k2) 5178 qi_mod_cas(l)= frac1*qi_prof_cas(k1) - frac2*qi_prof_cas(k2) 5179 u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2) 5180 v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2) 5181 ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2) 5182 vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2) 5183 w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2) 5184 omega_mod_cas(l)= frac1*omega_prof_cas(k1) - frac2*omega_prof_cas(k2) 5185 du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2) 5186 hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2) 5187 vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2) 5188 dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2) 5189 hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2) 5190 vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2) 5191 dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2) 5192 ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2) 5193 vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2) 5194 dth_mod_cas(l)= frac1*dth_prof_cas(k1) - frac2*dth_prof_cas(k2) 5195 hth_mod_cas(l)= frac1*hth_prof_cas(k1) - frac2*hth_prof_cas(k2) 5196 vth_mod_cas(l)= frac1*vth_prof_cas(k1) - frac2*vth_prof_cas(k2) 5197 dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2) 5198 hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2) 5199 vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2) 5200 5201 endif ! play.le.plev_prof_cas(1) 5202 5203 else ! above max altitude of forcing file 5204 5205 !jyg 5206 fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg 5207 fact = max(fact,0.) !jyg 5208 fact = exp(-fact) !jyg 5209 t_mod_cas(l)= t_prof_cas(nlev_cas) !jyg 5210 theta_mod_cas(l)= th_prof_cas(nlev_cas) !jyg 5211 thv_mod_cas(l)= thv_prof_cas(nlev_cas) !jyg 5212 thl_mod_cas(l)= thl_prof_cas(nlev_cas) !jyg 5213 qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact !jyg 5214 ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact !jyg 5215 qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact !jyg 5216 u_mod_cas(l)= u_prof_cas(nlev_cas)*fact !jyg 5217 v_mod_cas(l)= v_prof_cas(nlev_cas)*fact !jyg 5218 ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact !jyg 5219 vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact !jyg 5220 w_mod_cas(l)= 0.0 !jyg 5221 du_mod_cas(l)= du_prof_cas(nlev_cas)*fact 5222 hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact !jyg 5223 vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact !jyg 5224 dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact 5225 hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact !jyg 5226 vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact !jyg 5227 dt_mod_cas(l)= dt_prof_cas(nlev_cas) 5228 ht_mod_cas(l)= ht_prof_cas(nlev_cas) !jyg 5229 vt_mod_cas(l)= vt_prof_cas(nlev_cas) !jyg 5230 dth_mod_cas(l)= dth_prof_cas(nlev_cas) 5231 hth_mod_cas(l)= hth_prof_cas(nlev_cas) !jyg 5232 vth_mod_cas(l)= vth_prof_cas(nlev_cas) !jyg 5233 dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact 5234 hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact !jyg 5235 vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact !jyg 5236 5237 endif ! play 5238 5239 enddo ! l 5240 5241 ! do l = 1,llm 5242 ! print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ', 5243 ! $ l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) 5244 ! enddo 5245 5246 return 5247 end 5248 !***************************************************************************** 5249 5250 -
LMDZ5/branches/testing/libf/phylmd/dyn1d/1D_decl_cases.h
r2408 r2720 33 33 34 34 real w_mod(llm), t_mod(llm),q_mod(llm) 35 real u_mod(llm),v_mod(llm), ht_mod(llm),vt_mod(llm) 35 real u_mod(llm),v_mod(llm), ht_mod(llm),vt_mod(llm),ug_mod(llm),vg_mod(llm) 36 36 real hq_mod(llm),vq_mod(llm),qv_mod(llm),ql_mod(llm),qt_mod(llm) 37 37 real th_mod(llm) … … 94 94 95 95 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 96 !Declarations specifiques au cas GABLS4 (MPL 20141023) 97 character*80 :: fich_gabls4 98 integer nlev_gabls4, nt_gabls4, nsol_gabls4 99 parameter (nlev_gabls4=90, nt_gabls4=37, nsol_gabls4=19) 100 integer year_ini_gabls4, day_ini_gabls4, mth_ini_gabls4 101 real heure_ini_gabls4 102 real day_ju_ini_gabls4 ! Julian day of gabls4 first day 103 parameter (year_ini_gabls4=2009) 104 parameter (mth_ini_gabls4=12) 105 parameter (day_ini_gabls4=11) ! 11 = 11 decembre 2009 106 parameter (heure_ini_gabls4=0.) !0UTC en secondes 107 real dt_gabls4 108 parameter (dt_gabls4=3600.) ! 1 forcage ttes les heures 109 110 !profils initiaux: 111 real plev_gabls4(nlev_gabls4) 112 real zz_gabls4(nlev_gabls4) 113 real th_gabls4(nlev_gabls4),t_gabls4(nlev_gabls4),qv_gabls4(nlev_gabls4) 114 real u_gabls4(nlev_gabls4), v_gabls4(nlev_gabls4) 115 real depth_sn_gabls4(nsol_gabls4),tsnow_gabls4(nsol_gabls4),snow_dens_gabls4(nsol_gabls4) 116 real t_gabi(nlev_gabls4),qv_gabi(nlev_gabls4) 117 real u_gabi(nlev_gabls4), v_gabi(nlev_gabls4),ug_gabi(nlev_gabls4), vg_gabi(nlev_gabls4) 118 real ht_gabi(nlev_gabls4),hq_gabi(nlev_gabls4),poub(nlev_gabls4) 119 120 !forcings 121 real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4) 122 real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4) 123 real tg_gabls4(nt_gabls4) 124 real ht_profg(nlev_gabls4),hq_profg(nlev_gabls4) 125 real ug_profg(nlev_gabls4),vg_profg(nlev_gabls4) 126 real tg_profg 127 128 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 129 96 130 !Declarations specifiques au cas DICE (MPL 02072013) 97 131 character*80 :: fich_dice … … 112 146 113 147 real zz_dice(nlev_dice) 114 real t h_dice(nlev_dice),qv_dice(nlev_dice)148 real t_dice(nlev_dice),qv_dice(nlev_dice) 115 149 real u_dice(nlev_dice), v_dice(nlev_dice),o3_dice(nlev_dice) 116 150 real ht_dice(nlev_dice,nt_dice) … … 119 153 real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice) 120 154 real o3_mod(llm),hu_mod(llm),hv_mod(llm) 121 real t h_dicei(nlev_dice),qv_dicei(nlev_dice)155 real t_dicei(nlev_dice),qv_dicei(nlev_dice) 122 156 real u_dicei(nlev_dice), v_dicei(nlev_dice),o3_dicei(nlev_dice) 123 157 real ht_dicei(nlev_dice) … … 209 243 real thl_mod(llm),omega_mod(llm),o3mmr_mod(llm),tke_mod(llm) 210 244 !vertical advection computation 211 real d_t_z(llm), d_q_z(llm)212 real d_t_dyn_z(llm), d_q_dyn_z(llm)245 real d_t_z(llm),d_th_z(llm), d_q_z(llm) 246 real d_t_dyn_z(llm),d_th_dyn_z(llm), d_q_dyn_z(llm) 213 247 real d_u_z(llm),d_v_z(llm) 214 248 real d_u_dyn(llm),d_v_dyn(llm) … … 244 278 245 279 real w_mod_cas(llm), t_mod_cas(llm),q_mod_cas(llm) 280 real theta_mod_cas(llm),thl_mod_cas(llm),thv_mod_cas(llm) 281 real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm) 246 282 real ug_mod_cas(llm),vg_mod_cas(llm) 247 283 real u_mod_cas(llm),v_mod_cas(llm) 284 real omega_mod_cas(llm) 248 285 real ht_mod_cas(llm),vt_mod_cas(llm),dt_mod_cas(llm),dtrad_mod_cas(llm) 286 real hth_mod_cas(llm),vth_mod_cas(llm),dth_mod_cas(llm) 249 287 real hq_mod_cas(llm),vq_mod_cas(llm),dq_mod_cas(llm) 250 288 real hu_mod_cas(llm),vu_mod_cas(llm),du_mod_cas(llm) … … 253 291 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 254 292 293 -
LMDZ5/branches/testing/libf/phylmd/dyn1d/1D_interp_cases.h
r2594 r2720 118 118 ! vertical interpolation: 119 119 CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice & 120 & ,t h_dice,qv_dice,u_dice,v_dice,o3_dice &120 & ,t_dice,qv_dice,u_dice,v_dice,o3_dice & 121 121 & ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd,omega_profd & 122 & ,t h_mod,qv_mod,u_mod,v_mod,o3_mod &122 & ,t_mod,qv_mod,u_mod,v_mod,o3_mod & 123 123 & ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc) 124 124 ! do l = 1, llm … … 192 192 endif ! forcing_dice 193 193 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 194 ! Interpolation gabls4 forcing 195 !--------------------------------------------------------------------- 196 if (forcing_gabls4 ) then 197 198 if (prt_level.ge.1) then 199 print*,'#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_gabls4=',& 200 & day,day1,(day-day1)*86400.,(day-day1)*86400/dt_gabls4 201 endif 202 203 ! time interpolation: 204 CALL interp_gabls4_time(daytime,day1,annee_ref & 205 & ,year_ini_gabls4,day_ju_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4 & 206 & ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4 & 207 & ,ug_profg,vg_profg,ht_profg,hq_profg,tg_profg) 208 209 if (type_ts_forcing.eq.1) ts_cur = tg_prof ! SST used in read_tsurf1d 210 211 ! vertical interpolation: 212 ! on re-utilise le programme interp_dice_vertical: les transformations sur 213 ! plev_gabls4,th_gabls4,qv_gabls4,u_gabls4,v_gabls4 ne sont pas prises en compte. 214 ! seules celles sur ht_profg,hq_profg,ug_profg,vg_profg sont prises en compte. 215 216 CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4 & 217 ! & ,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,poub & 218 & ,poub,poub,poub,poub,poub & 219 & ,ht_profg,hq_profg,ug_profg,vg_profg,poub,poub & 220 & ,t_mod,qv_mod,u_mod,v_mod,o3_mod & 221 & ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc) 222 223 do l = 1, llm 224 ug(l)= ug_mod(l) 225 vg(l)= vg_mod(l) 226 d_th_adv(l)=ht_mod(l) 227 d_q_adv(l,1)=hq_mod(l) 228 enddo 229 230 endif ! forcing_gabls4 231 !--------------------------------------------------------------------- 232 194 233 !--------------------------------------------------------------------- 195 234 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 766 805 enddo 767 806 807 ! Faut-il multiplier par -1 ? (MPL 20160713) 808 IF(ok_flux_surf) THEN 809 fsens=sens_prof_cas 810 flat=lat_prof_cas 811 ENDIF 812 ! 813 IF (ok_prescr_ust) THEN 814 ust=ustar_prof_cas 815 print *,'ust=',ust 816 ENDIF 768 817 endif ! forcing_case 769 818 770 819 771 820 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 772 821 !--------------------------------------------------------------------- 822 ! Interpolation forcing standard case 823 !--------------------------------------------------------------------- 824 if (forcing_case2) then 825 826 print*, & 827 & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=', & 828 & daytime,day1,(daytime-day1)*86400., & 829 & (daytime-day1)*86400/pdt_cas 830 831 ! time interpolation: 832 CALL interp2_case_time(daytime,day1,annee_ref & 833 ! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas & 834 & ,nt_cas,nlev_cas & 835 & ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas & 836 & ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas & 837 & ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas & 838 & ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas & 839 & ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas & 840 ! 841 & ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas & 842 & ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas & 843 & ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas & 844 & ,du_prof_cas,hu_prof_cas,vu_prof_cas & 845 & ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas & 846 & ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & 847 & ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas & 848 & ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas) 849 850 ts_cur = ts_prof_cas 851 ! psurf=plev_prof_cas(1) 852 psurf=ps_prof_cas 853 854 ! vertical interpolation: 855 CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas & 856 & ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas & 857 & ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas & 858 & ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas & 859 & ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas & 860 & ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & 861 & ,dth_prof_cas,hth_prof_cas,vth_prof_cas & 862 ! 863 & ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas & 864 & ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas & 865 & ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas & 866 & ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas & 867 & ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc) 868 869 870 DO l=1,llm 871 teta(l)=temp(l)*(100000./play(l))**(rd/rcpd) 872 ENDDO 873 !calcul de l'advection verticale a partir du omega 874 !Calcul des gradients verticaux 875 !initialisation 876 d_t_z(:)=0. 877 d_th_z(:)=0. 878 d_q_z(:)=0. 879 d_t_dyn_z(:)=0. 880 d_th_dyn_z(:)=0. 881 d_q_dyn_z(:)=0. 882 DO l=2,llm-1 883 d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1)) 884 d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1)) 885 d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1)) 886 ENDDO 887 d_t_z(1)=d_t_z(2) 888 d_th_z(1)=d_th_z(2) 889 d_q_z(1)=d_q_z(2) 890 d_t_z(llm)=d_t_z(llm-1) 891 d_th_z(llm)=d_th_z(llm-1) 892 d_q_z(llm)=d_q_z(llm-1) 893 894 !Calcul de l advection verticale 895 d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:) 896 d_th_dyn_z(:)=w_mod_cas(:)*d_th_z(:) 897 d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:) 898 899 !wind nudging 900 if (nudging_u.gt.0.) then 901 do l=1,llm 902 u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u) 903 enddo 904 else 905 do l=1,llm 906 ug(l) = u_mod_cas(l) 907 enddo 908 endif 909 910 if (nudging_v.gt.0.) then 911 do l=1,llm 912 v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v) 913 enddo 914 else 915 do l=1,llm 916 vg(l) = v_mod_cas(l) 917 enddo 918 endif 919 920 if (nudging_w.gt.0.) then 921 do l=1,llm 922 w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w) 923 enddo 924 else 925 do l=1,llm 926 w(l) = w_mod_cas(l) 927 enddo 928 endif 929 930 !nudging of q and temp 931 if (nudging_t.gt.0.) then 932 do l=1,llm 933 temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t) 934 enddo 935 endif 936 if (nudging_q.gt.0.) then 937 do l=1,llm 938 q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q) 939 enddo 940 endif 941 942 do l = 1, llm 943 omega(l) = w_mod_cas(l) 944 omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 945 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 946 947 !calcul advection 948 ! if ((tend_u.eq.1).and.(tend_w.eq.0)) then 949 ! d_u_adv(l)=du_mod_cas(l) 950 ! else if ((tend_u.eq.1).and.(tend_w.eq.1)) then 951 ! d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l) 952 ! endif 953 ! 954 ! if ((tend_v.eq.1).and.(tend_w.eq.0)) then 955 ! d_v_adv(l)=dv_mod_cas(l) 956 ! else if ((tend_v.eq.1).and.(tend_w.eq.1)) then 957 ! d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l) 958 ! endif 959 ! 960 !----------------------------------------------------- 961 if (tadv.eq.1 .or. tadvh.eq.1) then 962 d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l) 963 else if (tadvv.eq.1) then 964 ! ATTENTION d_t_dyn_z pas calcule (voir twpice) 965 d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l) 966 endif 967 print *,'interp_case d_t_dyn_z=',d_t_dyn_z(l),d_q_dyn_z(l) 968 969 ! Verifier le signe !! 970 if (thadv.eq.1 .or. thadvh.eq.1) then 971 d_th_adv(l)=dth_mod_cas(l) 972 print *,'dthadv=',d_th_adv(l)*86400. 973 else if (thadvv.eq.1) then 974 d_th_adv(l)=hth_mod_cas(l)-d_th_dyn_z(l) 975 endif 976 977 ! Verifier le signe !! 978 if ((qadv.eq.1).and.(forc_w.eq.0)) then 979 d_q_adv(l,1)=dq_mod_cas(l) 980 else if ((qadvh.eq.1).and.(forc_w.eq.1)) then 981 d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l) 982 endif 983 984 if (trad.eq.1) then 985 tend_rayo=1 986 dt_cooling(l) = dtrad_mod_cas(l) 987 ! print *,'dt_cooling=',dt_cooling(l) 988 else 989 dt_cooling(l) = 0.0 990 endif 991 enddo 992 993 ! Faut-il multiplier par -1 ? (MPL 20160713) 994 IF(ok_flux_surf) THEN 995 fsens=-1.*sens_prof_cas 996 flat=-1.*lat_prof_cas 997 print *,'1D_interp: sens,flat',fsens,flat 998 ENDIF 999 ! 1000 IF (ok_prescr_ust) THEN 1001 ust=ustar_prof_cas 1002 print *,'ust=',ust 1003 ENDIF 1004 endif ! forcing_case2 1005 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1006 -
LMDZ5/branches/testing/libf/phylmd/dyn1d/1D_read_forc_cases.h
r2408 r2720 367 367 fich_dice='dice_driver.nc' 368 368 call read_dice(fich_dice,nlev_dice,nt_dice & 369 & ,zz_dice,plev_dice,t h_dice,qv_dice,u_dice,v_dice,o3_dice &369 & ,zz_dice,plev_dice,t_dice,qv_dice,u_dice,v_dice,o3_dice & 370 370 & ,shf_dice,lhf_dice,lwup_dice,swup_dice,tg_dice,ustar_dice& 371 371 & ,psurf_dice,ug_dice,vg_dice,ht_dice,hq_dice & … … 376 376 !champs initiaux: 377 377 do k=1,nlev_dice 378 t h_dicei(k)=th_dice(k)378 t_dicei(k)=t_dice(k) 379 379 qv_dicei(k)=qv_dice(k) 380 380 u_dicei(k)=u_dice(k) … … 405 405 406 406 CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice & 407 & ,t h_dicei,qv_dicei,u_dicei,v_dicei,o3_dicei &407 & ,t_dicei,qv_dicei,u_dicei,v_dicei,o3_dicei & 408 408 & ,ht_dicei,hq_dicei,hu_dicei,hv_dicei,w_dicei,omega_dicei& 409 & ,t h_mod,qv_mod,u_mod,v_mod,o3_mod &409 & ,t_mod,qv_mod,u_mod,v_mod,o3_mod & 410 410 & ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc) 411 411 … … 425 425 do l = 1, llm 426 426 ! Ligne du dessous ?? decommenter si on lit theta au lieu de temp 427 428 !temp(l) = t_mod(l)427 ! temp(l) = th_mod(l)*(play(l)/pzero)**rkappa 428 temp(l) = t_mod(l) 429 429 q(l,1) = qv_mod(l) 430 430 q(l,2) = 0.0 … … 473 473 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 474 474 !--------------------------------------------------------------------- 475 ! Forcing from GABLS4 experiment 476 !--------------------------------------------------------------------- 477 478 !!!! Si la temperature de surface n'est pas impos??e: 479 480 if (forcing_gabls4) then 481 !read GABLS4 forcings 482 483 fich_gabls4='gabls4_driver.nc' 484 485 486 call read_gabls4(fich_gabls4,nlev_gabls4,nt_gabls4,nsol_gabls4,zz_gabls4,depth_sn_gabls4,ug_gabls4,vg_gabls4 & 487 & ,plev_gabls4,th_gabls4,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,ht_gabls4,hq_gabls4,tg_gabls4,tsnow_gabls4,snow_dens_gabls4) 488 489 write(*,*) 'Forcing GABLS4 lu' 490 491 !champs initiaux: 492 do k=1,nlev_gabls4 493 t_gabi(k)=t_gabls4(k) 494 qv_gabi(k)=qv_gabls4(k) 495 u_gabi(k)=u_gabls4(k) 496 v_gabi(k)=v_gabls4(k) 497 poub(k)=0. 498 ht_gabi(k)=ht_gabls4(k,1) 499 hq_gabi(k)=hq_gabls4(k,1) 500 ug_gabi(k)=ug_gabls4(k,1) 501 vg_gabi(k)=vg_gabls4(k,1) 502 enddo 503 504 omega(:)=0. 505 omega2(:)=0. 506 rho(:)=0. 507 ! vertical interpolation using TOGA interpolation routine: 508 ! write(*,*)'avant interp vert', t_proftwp 509 ! 510 ! CALL interp_dice_time(daytime,day1,annee_ref 511 ! i ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice 512 ! i ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice 513 ! i ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice 514 ! i ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice 515 ! o ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof 516 ! o ,ustar_prof,psurf_prof,ug_profd,vg_profd 517 ! o ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd 518 ! o ,omega_profd) 519 520 CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4 & 521 & ,t_gabi,qv_gabi,u_gabi,v_gabi,poub & 522 & ,ht_gabi,hq_gabi,ug_gabi,vg_gabi,poub,poub & 523 & ,t_mod,qv_mod,u_mod,v_mod,o3_mod & 524 & ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc) 525 526 ! Les forcages GABLS4 ont l air d etre en K/S quoiqu en dise le fichier gabls4_driver.nc !? MPL 20141024 527 ! ht_mod(:)=ht_mod(:)/86400. 528 ! hq_mod(:)=hq_mod(:)/86400. 529 530 ! initial and boundary conditions : 531 write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc 532 do l = 1, llm 533 ! Ligne du dessous ?? decommenter si on lit theta au lieu de temp 534 ! temp(l) = th_mod(l)*(play(l)/pzero)**rkappa 535 temp(l) = t_mod(l) 536 q(l,1) = qv_mod(l) 537 q(l,2) = 0.0 538 ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) 539 u(l) = u_mod(l) 540 v(l) = v_mod(l) 541 ug(l)=ug_mod(l) 542 vg(l)=vg_mod(l) 543 544 ! 545 ! tg=tsurf 546 ! 547 548 print *,'***** tsurf=',tsurf 549 rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 550 ! omega(l) = w_mod(l)*(-rg*rho(l)) 551 omega(l) = omega_mod(l) 552 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 553 554 555 556 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 557 !on applique le forcage total au premier pas de temps 558 !attention: signe different de toga 559 ! d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) 560 !forcage en th 561 d_th_adv(l) = ht_mod(l) 562 d_q_adv(l,1) = hq_mod(l) 563 d_q_adv(l,2) = 0.0 564 dt_cooling(l)=0. 565 enddo 566 567 !--------------- Residus forcages du cas Dice (a supprimer) MPL 20141024--------------- 568 ! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par 569 ! le coefficient de trainee en surface cd**2=ustar*vent(k=1) 570 ! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface 571 ! MPL 05082013 572 ! ust=ustar_dice(1) 573 ! tg=tg_dice(1) 574 ! print *,'ust= ',ust 575 ! IF (tsurf .LE. 0.) THEN 576 ! tsurf= tg_dice(1) 577 ! ENDIF 578 ! psurf= psurf_dice(1) 579 ! solsw_in = (1.-albedo)/albedo*swup_dice(1) 580 ! sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1) 581 ! PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in 582 !-------------------------------------------------------------------------------------- 583 endif !forcing_gabls4 584 585 586 475 587 ! Forcing from Arm_Cu case 476 588 ! For this case, ifa_armcu.txt contains sensible, latent heat fluxes … … 797 909 endif !forcing_case 798 910 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 911 !--------------------------------------------------------------------- 912 ! Forcing from standard case : 913 !--------------------------------------------------------------------- 914 915 if (forcing_case2) then 916 917 write(*,*),'avant call read2_1D_cas' 918 call read2_1D_cas 919 write(*,*) 'Forcing read' 920 921 !Time interpolation for initial conditions using interpolation routine 922 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1 923 CALL interp2_case_time(daytime,day1,annee_ref & 924 ! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas & 925 & ,nt_cas,nlev_cas & 926 & ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas & 927 & ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas & 928 & ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas & 929 & ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas & 930 & ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas & 931 ! 932 & ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas & 933 & ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas & 934 & ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas & 935 & ,du_prof_cas,hu_prof_cas,vu_prof_cas & 936 & ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas & 937 & ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & 938 & ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas & 939 & ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas) 940 941 do l = 1, nlev_cas 942 print *,'apres 1ere interp: plev_cas, plev_prof_cas=',l,plev_cas(l,1),plev_prof_cas(l) 943 enddo 944 945 ! vertical interpolation using interpolation routine: 946 ! write(*,*)'avant interp vert', t_prof 947 CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas & 948 & ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas & 949 & ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas & 950 & ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas & 951 & ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas & 952 & ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & 953 & ,dth_prof_cas,hth_prof_cas,vth_prof_cas & 954 ! 955 & ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas & 956 & ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas & 957 & ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas & 958 & ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas & 959 & ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc) 960 961 ! write(*,*) 'Profil initial forcing case interpole',t_mod 962 963 ! initial and boundary conditions : 964 ! tsurf = ts_prof_cas 965 ts_cur = ts_prof_cas 966 psurf=plev_prof_cas(1) 967 write(*,*) 'SST initiale: ',tsurf 968 do l = 1, llm 969 temp(l) = t_mod_cas(l) 970 q(l,1) = qv_mod_cas(l) 971 q(l,2) = ql_mod_cas(l) 972 u(l) = u_mod_cas(l) 973 ug(l)= u_mod_cas(l) 974 v(l) = v_mod_cas(l) 975 vg(l)= v_mod_cas(l) 976 omega(l) = w_mod_cas(l) 977 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 978 979 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 980 !on applique le forcage total au premier pas de temps 981 !attention: signe different de toga 982 d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l)) 983 d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l)) 984 ! d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l)) 985 d_q_adv(l,1) = dq_mod_cas(l) 986 d_q_adv(l,2) = 0.0 987 ! d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l)) 988 d_u_adv(l) = du_mod_cas(l) 989 ! d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l)) 990 d_u_adv(l) = dv_mod_cas(l) 991 enddo 992 993 ! Faut-il multiplier par -1 ? (MPL 20160713) 994 IF (ok_flux_surf) THEN 995 fsens=-1.*sens_prof_cas 996 flat=-1.*lat_prof_cas 997 ENDIF 998 ! 999 IF (ok_prescr_ust) THEN 1000 ust=ustar_prof_cas 1001 print *,'ust=',ust 1002 ENDIF 1003 1004 endif !forcing_case2 1005 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1006 -
LMDZ5/branches/testing/libf/phylmd/dyn1d/compar1d.h
r2408 r2720 9 9 real :: tsurf 10 10 real :: rugos 11 real :: rugosh 11 12 real :: xqsol(1:2) 12 13 real :: qsurf … … 14 15 real :: zsurf 15 16 real :: albedo 17 real :: snowmass 16 18 17 19 real :: time … … 30 32 logical :: ok_old_disvert 31 33 34 ! Pour les forcages communs: ces entiers valent 0 ou 1 35 ! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale 36 ! idem pour l advection en theta 37 ! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale 38 ! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv) 39 ! forcages en omega, w, vent geostrophique ou ustar 40 ! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging 41 42 integer :: tadv, tadvv, tadvh, qadv, qadvv, qadvh, thadv, thadvv, thadvh, trad 43 integer :: forc_omega, forc_w, forc_geo, forc_ustar 44 real :: nudging_u, nudging_v, nudging_w, nudging_t, nudging_q 32 45 common/com_par1d/ & 33 & nat_surf,tsurf,rugos, 46 & nat_surf,tsurf,rugos,rugosh, & 34 47 & xqsol,qsurf,psurf,zsurf,albedo,time,time_ini,xlat,xlon,airefi, & 35 48 & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp, & 36 49 & forcing_type,tend_u,tend_v,tend_w,tend_t,tend_q,tend_rayo, & 37 50 & nudge_u,nudge_v,nudge_w,nudge_t,nudge_q, & 38 & iflag_nudge, & 39 & restart,ok_old_disvert 51 & iflag_nudge,snowmass, & 52 & restart,ok_old_disvert, & 53 & tadv, tadvv, tadvh, qadv, qadvv, qadvh, thadv, thadvv, thadvh, & 54 & trad, forc_omega, forc_w, forc_geo, forc_ustar, & 55 & nudging_u, nudging_v, nudging_t, nudging_q 40 56 41 57 !$OMP THREADPRIVATE(/com_par1d/) … … 50 66 51 67 68 -
LMDZ5/branches/testing/libf/phylmd/dyn1d/lmdz1d.F90
r2641 r2720 21 21 zgam, zmax0, zmea, zpic, zsig, & 22 22 zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl 23 23 24 USE dimphy 24 25 USE surface_data, only : type_ocean,ok_veget … … 31 32 USE indice_sol_mod 32 33 USE phyaqua_mod 33 USE mod_1D_cases_read 34 ! USE mod_1D_cases_read 35 USE mod_1D_cases_read2 34 36 USE mod_1D_amma_read 35 37 USE print_control_mod, ONLY: lunout, prt_level … … 131 133 logical :: forcing_amma = .false. 132 134 logical :: forcing_dice = .false. 135 logical :: forcing_gabls4 = .false. 136 133 137 logical :: forcing_GCM2SCM = .false. 134 138 logical :: forcing_GCSSold = .false. … … 137 141 logical :: forcing_fire = .false. 138 142 logical :: forcing_case = .false. 143 logical :: forcing_case2 = .false. 139 144 integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file 140 145 ! (cf read_tsurf1d.F) … … 174 179 real :: pzero=1.e5 175 180 real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1) 176 real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1) ,poub181 real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1) 177 182 178 183 !--------------------------------------------------------------------- … … 189 194 real :: du_phys(llm),dv_phys(llm),dt_phys(llm) 190 195 real :: dt_dyn(llm) 191 real :: dt_cooling(llm),d_t h_adv(llm),d_t_nudge(llm)196 real :: dt_cooling(llm),d_t_adv(llm),d_th_adv(llm),d_t_nudge(llm) 192 197 real :: d_u_nudge(llm),d_v_nudge(llm) 193 198 real :: du_adv(llm),dv_adv(llm) … … 322 327 ! Different stages: soil model alone, atm. model alone 323 328 ! then both models coupled 329 !forcing_type = 8 ==> forcing_gabls4 = .true. 330 ! initial profiles and large scale forcings in gabls4_driver.nc 324 331 !forcing_type >= 100 ==> forcing_case = .true. 325 332 ! initial profiles and large scale forcings in cas.nc … … 327 334 ! 101=cindynamo 328 335 ! 102=bomex 336 !forcing_type >= 100 ==> forcing_case2 = .true. 337 ! temporary flag while all the 1D cases are not whith the same cas.nc forcing file 338 ! 103=arm_cu2 ie arm_cu with new forcing format 339 ! 104=rico2 ie rico with new forcing format 329 340 !forcing_type = 40 ==> forcing_GCSSold = .true. 330 341 ! initial profile from GCSS file … … 363 374 elseif (forcing_type .eq.7) THEN 364 375 forcing_dice = .true. 376 elseif (forcing_type .eq.8) THEN 377 forcing_gabls4 = .true. 365 378 elseif (forcing_type .eq.101) THEN ! Cindynamo starts 1-10-2011 0h 366 379 forcing_case = .true. … … 375 388 mth_ini_cas=6 376 389 day_deb=24 390 heure_ini_cas=0. 391 pdt_cas=1800. ! forcing frequency 392 elseif (forcing_type .eq.103) THEN ! Arm_cu starts 21-6-1997 11h30 393 forcing_case2 = .true. 394 year_ini_cas=1997 395 mth_ini_cas=6 396 day_deb=21 397 heure_ini_cas=11.5 398 pdt_cas=1800. ! forcing frequency 399 elseif (forcing_type .eq.104) THEN ! rico starts 16-12-2004 0h 400 forcing_case2 = .true. 401 year_ini_cas=2004 402 mth_ini_cas=12 403 day_deb=16 377 404 heure_ini_cas=0. 378 405 pdt_cas=1800. ! forcing frequency … … 449 476 endif 450 477 print *,'fnday=',fnday 451 478 ! start_time doit etre en FRACTION DE JOUR 452 479 start_time=time_ini/24. 453 480 454 481 ! Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026) 455 482 IF(forcing_type .EQ. 61) fnday=53100./86400. 483 IF(forcing_type .EQ. 103) fnday=53100./86400. 456 484 ! Special case for amma which lasts less than one day : 64800s !! (MPL 20120216) 457 485 IF(forcing_type .EQ. 6) fnday=64800./86400. 458 486 ! IF(forcing_type .EQ. 6) fnday=50400./86400. 487 IF(forcing_type .EQ. 8 ) fnday=129600./86400. 459 488 annee_ref = anneeref 460 489 mois = 1 … … 487 516 & (year_ini_dice,mth_ini_dice,day_ini_dice,heure_ini_dice & 488 517 & ,day_ju_ini_dice) 518 ELSEIF (forcing_type .eq.8 ) THEN 519 ! Convert the initial date of GABLS4 to Julian day 520 call ymds2ju & 521 & (year_ini_gabls4,mth_ini_gabls4,day_ini_gabls4,heure_ini_gabls4 & 522 & ,day_ju_ini_gabls4) 489 523 ELSEIF (forcing_type .gt.100) THEN 490 524 ! Convert the initial date to Julian day … … 492 526 print*,'time case',year_ini_cas,mth_ini_cas,day_ini_cas 493 527 call ymds2ju & 494 & (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas 528 & (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas*3600 & 495 529 & ,day_ju_ini_cas) 496 530 print*,'time case 2',day_ini_cas,day_ju_ini_cas … … 514 548 ENDIF 515 549 550 IF (forcing_type .gt.100) THEN 551 daytime = day + heure_ini_cas/24. ! 1st day and initial time of the simulation 552 ELSE 516 553 daytime = day + time_ini/24. ! 1st day and initial time of the simulation 554 ENDIF 517 555 ! Print out the actual date of the beginning of the simulation : 518 556 call ju2ymds(daytime,year_print, month_print,day_print,sec_print) … … 699 737 700 738 fder=0. 701 snsrf(1,:)= 0. ! couverture de neige des sous surface739 snsrf(1,:)=snowmass ! masse de neige des sous surface 702 740 qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface 703 741 fevap=0. 704 742 z0m(1,:)=rugos ! couverture de neige des sous surface 705 z0h(1,:)=rugos 743 z0h(1,:)=rugosh ! couverture de neige des sous surface 706 744 agesno = xagesno 707 745 tsoil(:,:,:)=tsurf … … 726 764 print*,'avant phyredem' 727 765 pctsrf(1,:)=0. 728 if (nat_surf.eq.0.) then766 if (nat_surf.eq.0.) then 729 767 pctsrf(1,is_oce)=1. 730 768 pctsrf(1,is_ter)=0. 731 else 769 pctsrf(1,is_lic)=0. 770 pctsrf(1,is_sic)=0. 771 else if (nat_surf .eq. 1) then 732 772 pctsrf(1,is_oce)=0. 733 773 pctsrf(1,is_ter)=1. 734 end if 774 pctsrf(1,is_lic)=0. 775 pctsrf(1,is_sic)=0. 776 else if (nat_surf .eq. 2) then 777 pctsrf(1,is_oce)=0. 778 pctsrf(1,is_ter)=0. 779 pctsrf(1,is_lic)=1. 780 pctsrf(1,is_sic)=0. 781 else if (nat_surf .eq. 3) then 782 pctsrf(1,is_oce)=0. 783 pctsrf(1,is_ter)=0. 784 pctsrf(1,is_lic)=0. 785 pctsrf(1,is_sic)=1. 786 787 end if 788 735 789 736 790 print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf & … … 1005 1059 1006 1060 if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice & 1007 & .or.forcing_amma ) then1061 & .or.forcing_amma .or. forcing_type.eq.101) then 1008 1062 fcoriolis=0.0 ; ug=0. ; vg=0. 1009 1063 endif 1010 if(forcing_rico) then 1064 1065 if(forcing_rico) then 1011 1066 dt_cooling=0. 1012 1067 endif 1013 1068 1014 1069 IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', & … … 1172 1227 !#endif 1173 1228 1229
Note: See TracChangeset
for help on using the changeset viewer.
