Ignore:
Timestamp:
Oct 17, 2016, 9:47:47 AM (8 years ago)
Author:
fhourdin
Message:

Introduction du cas GABLS4 par Etienne Vignon

Location:
LMDZ5/trunk/libf/phylmd/dyn1d
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/dyn1d/1DUTILS.h

    r2668 r2672  
    303303       CALL getin('rugos',rugos)
    304304
     305!Config  Key  = rugosh
     306!Config  Desc = coefficient de frottement
     307!Config  Def  = rugos
     308!Config  Help = calcul du Cdrag
     309       rugosh = rugos
     310       CALL getin('rugosh',rugosh)
     311
     312
     313
     314!Config  Key  = snowmass
     315!Config  Desc = mass de neige de la surface en kg/m2
     316!Config  Def  = 0.0000
     317!Config  Help = snowmass
     318       snowmass = 0.0000
     319       CALL getin('snowmass',snowmass)
     320
    305321!Config  Key  = wtsurf et wqsurf
    306322!Config  Desc = ???
    … …  
    395411      write(lunout,*)' zsurf = ', zsurf
    396412      write(lunout,*)' rugos = ', rugos
     413      write(lunout,*)' snowmass=', snowmass
    397414      write(lunout,*)' wtsurf = ', wtsurf
    398415      write(lunout,*)' wqsurf = ', wqsurf
    … …  
    31413158
    31423159!======================================================================
     3160        SUBROUTINE interp_gabls4_time(day,day1,annee_ref                              &
     3161     &             ,year_ini_gabls4,day_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4    &
     3162     &             ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4                          &
     3163     &             ,ug_prof,vg_prof,ht_prof,hq_prof,tg_prof)
     3164        implicit none
     3165
     3166!---------------------------------------------------------------------------------------
     3167! Time interpolation of a 2D field to the timestep corresponding to day
     3168!
     3169! day: current julian day
     3170! day1: first day of the simulation
     3171! nt_gabls4: total nb of data in the forcing (e.g. 37 for gabls4)
     3172! dt_gabls4: total time interval (in sec) between 2 forcing data (e.g. 60min. for gabls4)
     3173!---------------------------------------------------------------------------------------
     3174
     3175#include "compar1d.h"
     3176
     3177! inputs:
     3178        integer annee_ref
     3179        integer nt_gabls4,nlev_gabls4
     3180        integer year_ini_gabls4
     3181        real day, day1,day_ini_gabls4,dt_gabls4
     3182        real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4)
     3183        real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4)
     3184        real tg_gabls4(nt_gabls4), tg_prof
     3185! outputs:
     3186        real ug_prof(nlev_gabls4),vg_prof(nlev_gabls4)
     3187        real ht_prof(nlev_gabls4),hq_prof(nlev_gabls4)
     3188! local:
     3189        integer it_gabls41, it_gabls42,k
     3190        real timeit,time_gabls41,time_gabls42,frac
     3191
     3192
     3193
     3194! Check that initial day of the simulation consistent with gabls4 period:
     3195       if (forcing_type.eq.8 ) then
     3196       print *,'annee_ref=',annee_ref
     3197       print *,'day1=',day1
     3198       print *,'day_ini_gabls4=',day_ini_gabls4
     3199       if (annee_ref.ne.2009) then
     3200        print*,'Pour gabls4, annee_ref doit etre 2009'
     3201        print*,'Changer annee_ref dans run.def'
     3202        stop
     3203       endif
     3204       if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4) then
     3205        print*,'gabls4 a debute le 11 dec 2009 (jour julien=345)'
     3206        print*,'Changer dayref dans run.def',day1,day_ini_gabls4
     3207        stop
     3208       endif
     3209       if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4+2) then
     3210        print*,'gabls4 a fini le 12 dec 2009 (jour julien=346)'
     3211        print*,'Changer dayref ou nday dans run.def',day1,day_ini_gabls4
     3212        stop
     3213       endif
     3214       endif
     3215
     3216      timeit=(day-day_ini_gabls4)*86400
     3217       print *,'day,day_ini_gabls4=',day,day_ini_gabls4
     3218       print *,'nt_gabls4,dt,timeit=',nt_gabls4,dt_gabls4,timeit
     3219
     3220! Determine the closest observation times:
     3221       it_gabls41=INT(timeit/dt_gabls4)+1
     3222       it_gabls42=it_gabls41 + 1
     3223       time_gabls41=(it_gabls41-1)*dt_gabls4
     3224       time_gabls42=(it_gabls42-1)*dt_gabls4
     3225
     3226       if (it_gabls41 .ge. nt_gabls4) then
     3227        write(*,*) 'PB-stop: day, it_gabls41, it_gabls42, timeit: ',day,it_gabls41,it_gabls42,timeit/86400.
     3228        stop
     3229       endif
     3230
     3231! time interpolation:
     3232       frac=(time_gabls42-timeit)/(time_gabls42-time_gabls41)
     3233       frac=max(frac,0.0)
     3234
     3235
     3236       do k=1,nlev_gabls4
     3237        ug_prof(k) = ug_gabls4(k,it_gabls42)-frac*(ug_gabls4(k,it_gabls42)-ug_gabls4(k,it_gabls41))
     3238        vg_prof(k) = vg_gabls4(k,it_gabls42)-frac*(vg_gabls4(k,it_gabls42)-vg_gabls4(k,it_gabls41))
     3239        ht_prof(k) = ht_gabls4(k,it_gabls42)-frac*(ht_gabls4(k,it_gabls42)-ht_gabls4(k,it_gabls41))
     3240        hq_prof(k) = hq_gabls4(k,it_gabls42)-frac*(hq_gabls4(k,it_gabls42)-hq_gabls4(k,it_gabls41))
     3241        enddo
     3242        tg_prof=tg_gabls4(it_gabls42)-frac*(tg_gabls4(it_gabls42)-tg_gabls4(it_gabls41))
     3243        return
     3244        END
     3245
     3246!======================================================================
    31433247        SUBROUTINE interp_armcu_time(day,day1,annee_ref                    &
    31443248     &             ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu         &
    … …  
    40954199         end subroutine read_dice
    40964200!=====================================================================
     4201      subroutine read_gabls4(fich_gabls4,nlevel,ntime,nsol                    &
     4202     &     ,zz,depth_sn,ug,vg,pf,th,t,qv,u,v,hadvt,hadvq,tg,tsnow,snow_dens)
     4203
     4204!program reading initial profils and forcings of the Gabls4 case study
     4205
     4206
     4207      implicit none
     4208
     4209#include "netcdf.inc"
     4210
     4211      integer ntime,nlevel,nsol
     4212      integer l,k
     4213      character*80 :: fich_gabls4
     4214      real*8 time(ntime)
     4215
     4216!  ATTENTION: visiblement quand on lit gabls4_driver.nc on recupere les donnees
     4217! dans un ordre inverse par rapport a la convention LMDZ
     4218! ==> il faut tout inverser  (MPL 20141024)
     4219! les variables indexees "_i" sont celles qui sont lues dans gabls4_driver.nc
     4220      real*8 zz_i(nlevel),th_i(nlevel),pf_i(nlevel),t_i(nlevel)
     4221      real*8 qv_i(nlevel),u_i(nlevel),v_i(nlevel),ug_i(nlevel,ntime),vg_i(nlevel,ntime)
     4222      real*8 hadvt_i(nlevel,ntime),hadvq_i(nlevel,ntime)
     4223
     4224      real*8 zz(nlevel),th(nlevel),pf(nlevel),t(nlevel)
     4225      real*8 qv(nlevel),u(nlevel),v(nlevel),ug(nlevel,ntime),vg(nlevel,ntime)
     4226      real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime)
     4227
     4228      real*8 depth_sn(nsol),tsnow(nsol),snow_dens(nsol)
     4229      real*8 tg(ntime)
     4230      integer nid, ierr
     4231      integer nbvar3d
     4232      parameter(nbvar3d=30)
     4233      integer var3didin(nbvar3d)
     4234
     4235      ierr = NF_OPEN(fich_gabls4,NF_NOWRITE,nid)
     4236      if (ierr.NE.NF_NOERR) then
     4237         write(*,*) 'ERROR: Pb opening forcings nc file '
     4238         write(*,*) NF_STRERROR(ierr)
     4239         stop ""
     4240      endif
     4241
     4242
     4243       ierr=NF_INQ_VARID(nid,"height",var3didin(1))
     4244         if(ierr/=NF_NOERR) then
     4245           write(*,*) NF_STRERROR(ierr)
     4246           stop 'height'
     4247         endif
     4248
     4249      ierr=NF_INQ_VARID(nid,"depth_sn",var3didin(2))
     4250         if(ierr/=NF_NOERR) then
     4251           write(*,*) NF_STRERROR(ierr)
     4252           stop 'depth_sn'
     4253      endif
     4254
     4255      ierr=NF_INQ_VARID(nid,"Ug",var3didin(3))
     4256         if(ierr/=NF_NOERR) then
     4257           write(*,*) NF_STRERROR(ierr)
     4258           stop 'Ug'
     4259      endif
     4260
     4261      ierr=NF_INQ_VARID(nid,"Vg",var3didin(4))
     4262         if(ierr/=NF_NOERR) then
     4263           write(*,*) NF_STRERROR(ierr)
     4264           stop 'Vg'
     4265      endif
     4266       ierr=NF_INQ_VARID(nid,"pf",var3didin(5))
     4267         if(ierr/=NF_NOERR) then
     4268           write(*,*) NF_STRERROR(ierr)
     4269           stop 'pf'
     4270         endif
     4271
     4272      ierr=NF_INQ_VARID(nid,"theta",var3didin(6))
     4273         if(ierr/=NF_NOERR) then
     4274           write(*,*) NF_STRERROR(ierr)
     4275           stop 'theta'
     4276         endif
     4277
     4278      ierr=NF_INQ_VARID(nid,"tempe",var3didin(7))
     4279         if(ierr/=NF_NOERR) then
     4280           write(*,*) NF_STRERROR(ierr)
     4281           stop 'tempe'
     4282         endif
     4283
     4284      ierr=NF_INQ_VARID(nid,"qv",var3didin(8))
     4285         if(ierr/=NF_NOERR) then
     4286           write(*,*) NF_STRERROR(ierr)
     4287           stop 'qv'
     4288         endif
     4289
     4290      ierr=NF_INQ_VARID(nid,"u",var3didin(9))
     4291         if(ierr/=NF_NOERR) then
     4292           write(*,*) NF_STRERROR(ierr)
     4293           stop 'u'
     4294         endif
     4295
     4296      ierr=NF_INQ_VARID(nid,"v",var3didin(10))
     4297         if(ierr/=NF_NOERR) then
     4298           write(*,*) NF_STRERROR(ierr)
     4299           stop 'v'
     4300         endif
     4301
     4302      ierr=NF_INQ_VARID(nid,"hadvT",var3didin(11))
     4303         if(ierr/=NF_NOERR) then
     4304           write(*,*) NF_STRERROR(ierr)
     4305           stop 'hadvt'
     4306         endif
     4307
     4308      ierr=NF_INQ_VARID(nid,"hadvQ",var3didin(12))
     4309         if(ierr/=NF_NOERR) then
     4310           write(*,*) NF_STRERROR(ierr)
     4311           stop 'hadvq'
     4312      endif
     4313
     4314      ierr=NF_INQ_VARID(nid,"Tsnow",var3didin(14))
     4315         if(ierr/=NF_NOERR) then
     4316           write(*,*) NF_STRERROR(ierr)
     4317           stop 'tsnow'
     4318      endif
     4319
     4320      ierr=NF_INQ_VARID(nid,"snow_density",var3didin(15))
     4321         if(ierr/=NF_NOERR) then
     4322           write(*,*) NF_STRERROR(ierr)
     4323           stop 'snow_density'
     4324      endif
     4325
     4326      ierr=NF_INQ_VARID(nid,"Tg",var3didin(16))
     4327         if(ierr/=NF_NOERR) then
     4328           write(*,*) NF_STRERROR(ierr)
     4329           stop 'Tg'
     4330      endif
     4331
     4332
     4333!dimensions lecture
     4334!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
     4335 
     4336#ifdef NC_DOUBLE
     4337         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz_i)
     4338#else
     4339         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz_i)
     4340#endif
     4341         if(ierr/=NF_NOERR) then
     4342            write(*,*) NF_STRERROR(ierr)
     4343            stop "getvarup"
     4344         endif
     4345 
     4346#ifdef NC_DOUBLE
     4347         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),depth_sn)
     4348#else
     4349         ierr = NF_GET_VAR_REAL(nid,var3didin(2),depth_sn)
     4350#endif
     4351         if(ierr/=NF_NOERR) then
     4352            write(*,*) NF_STRERROR(ierr)
     4353            stop "getvarup"
     4354         endif
     4355 
     4356#ifdef NC_DOUBLE
     4357         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),ug_i)
     4358#else
     4359         ierr = NF_GET_VAR_REAL(nid,var3didin(3),ug_i)
     4360#endif
     4361         if(ierr/=NF_NOERR) then
     4362            write(*,*) NF_STRERROR(ierr)
     4363            stop "getvarup"
     4364         endif
     4365 
     4366#ifdef NC_DOUBLE
     4367         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),vg_i)
     4368#else
     4369         ierr = NF_GET_VAR_REAL(nid,var3didin(4),vg_i)
     4370#endif
     4371         if(ierr/=NF_NOERR) then
     4372            write(*,*) NF_STRERROR(ierr)
     4373            stop "getvarup"
     4374         endif
     4375 
     4376#ifdef NC_DOUBLE
     4377         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),pf_i)
     4378#else
     4379         ierr = NF_GET_VAR_REAL(nid,var3didin(5),pf_i)
     4380#endif
     4381         if(ierr/=NF_NOERR) then
     4382            write(*,*) NF_STRERROR(ierr)
     4383            stop "getvarup"
     4384         endif
     4385
     4386#ifdef NC_DOUBLE
     4387         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),th_i)
     4388#else
     4389         ierr = NF_GET_VAR_REAL(nid,var3didin(6),th_i)
     4390#endif
     4391         if(ierr/=NF_NOERR) then
     4392            write(*,*) NF_STRERROR(ierr)
     4393            stop "getvarup"
     4394         endif
     4395
     4396#ifdef NC_DOUBLE
     4397         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),t_i)
     4398#else
     4399         ierr = NF_GET_VAR_REAL(nid,var3didin(7),t_i)
     4400#endif
     4401         if(ierr/=NF_NOERR) then
     4402            write(*,*) NF_STRERROR(ierr)
     4403            stop "getvarup"
     4404         endif
     4405
     4406#ifdef NC_DOUBLE
     4407         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),qv_i)
     4408#else
     4409         ierr = NF_GET_VAR_REAL(nid,var3didin(8),qv_i)
     4410#endif
     4411         if(ierr/=NF_NOERR) then
     4412            write(*,*) NF_STRERROR(ierr)
     4413            stop "getvarup"
     4414         endif
     4415 
     4416#ifdef NC_DOUBLE
     4417         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),u_i)
     4418#else
     4419         ierr = NF_GET_VAR_REAL(nid,var3didin(9),u_i)
     4420#endif
     4421         if(ierr/=NF_NOERR) then
     4422            write(*,*) NF_STRERROR(ierr)
     4423            stop "getvarup"
     4424         endif
     4425 
     4426#ifdef NC_DOUBLE
     4427         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),v_i)
     4428#else
     4429         ierr = NF_GET_VAR_REAL(nid,var3didin(10),v_i)
     4430#endif
     4431         if(ierr/=NF_NOERR) then
     4432            write(*,*) NF_STRERROR(ierr)
     4433            stop "getvarup"
     4434         endif
     4435 
     4436#ifdef NC_DOUBLE
     4437         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),hadvt_i)
     4438#else
     4439         ierr = NF_GET_VAR_REAL(nid,var3didin(11),hadvt_i)
     4440#endif
     4441         if(ierr/=NF_NOERR) then
     4442            write(*,*) NF_STRERROR(ierr)
     4443            stop "getvarup"
     4444         endif
     4445 
     4446#ifdef NC_DOUBLE
     4447         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),hadvq_i)
     4448#else
     4449         ierr = NF_GET_VAR_REAL(nid,var3didin(12),hadvq_i)
     4450#endif
     4451         if(ierr/=NF_NOERR) then
     4452            write(*,*) NF_STRERROR(ierr)
     4453            stop "getvarup"
     4454         endif
     4455 
     4456#ifdef NC_DOUBLE
     4457         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),tsnow)
     4458#else
     4459         ierr = NF_GET_VAR_REAL(nid,var3didin(14),tsnow)
     4460#endif
     4461         if(ierr/=NF_NOERR) then
     4462            write(*,*) NF_STRERROR(ierr)
     4463            stop "getvarup"
     4464         endif
     4465 
     4466#ifdef NC_DOUBLE
     4467         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),snow_dens)
     4468#else
     4469         ierr = NF_GET_VAR_REAL(nid,var3didin(15),snow_dens)
     4470#endif
     4471         if(ierr/=NF_NOERR) then
     4472            write(*,*) NF_STRERROR(ierr)
     4473            stop "getvarup"
     4474         endif
     4475
     4476#ifdef NC_DOUBLE
     4477         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),tg)
     4478#else
     4479         ierr = NF_GET_VAR_REAL(nid,var3didin(16),tg)
     4480#endif
     4481         if(ierr/=NF_NOERR) then
     4482            write(*,*) NF_STRERROR(ierr)
     4483            stop "getvarup"
     4484         endif
     4485
     4486! On remet les variables lues dans le bon ordre des niveaux (MPL 20141024)
     4487         do k=1,nlevel
     4488           zz(k)=zz_i(nlevel+1-k)
     4489           ug(k,:)=ug_i(nlevel+1-k,:)
     4490           vg(k,:)=vg_i(nlevel+1-k,:)
     4491           pf(k)=pf_i(nlevel+1-k)
     4492           print *,'pf=',pf(k)
     4493           th(k)=th_i(nlevel+1-k)
     4494           t(k)=t_i(nlevel+1-k)
     4495           qv(k)=qv_i(nlevel+1-k)
     4496           u(k)=u_i(nlevel+1-k)
     4497           v(k)=v_i(nlevel+1-k)
     4498           hadvt(k,:)=hadvt_i(nlevel+1-k,:)
     4499           hadvq(k,:)=hadvq_i(nlevel+1-k,:)
     4500         enddo
     4501         return
     4502 end subroutine read_gabls4
     4503!=====================================================================
     4504
    40974505!     Reads CIRC input files     
    40984506

  • LMDZ5/trunk/libf/phylmd/dyn1d/1D_decl_cases.h

    r2373 r2672  
    3333
    3434        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)
    3636        real hq_mod(llm),vq_mod(llm),qv_mod(llm),ql_mod(llm),qt_mod(llm)
    3737        real th_mod(llm)
    … …  
    9494
    9595!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     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
    96130!Declarations specifiques au cas DICE     (MPL 02072013)
    97131        character*80 :: fich_dice

  • LMDZ5/trunk/libf/phylmd/dyn1d/1D_interp_cases.h

    r2565 r2672  
    192192      endif ! forcing_dice
    193193!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     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
    194233!---------------------------------------------------------------------
    195234!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  • LMDZ5/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h

    r2332 r2672  
    473473!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    474474!---------------------------------------------------------------------
     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
    475587! Forcing from Arm_Cu case                   
    476588! For this case, ifa_armcu.txt contains sensible, latent heat fluxes

  • LMDZ5/trunk/libf/phylmd/dyn1d/compar1d.h

    r2328 r2672  
    99      real :: tsurf
    1010      real :: rugos
     11      real :: rugosh
    1112      real :: xqsol(1:2)
    1213      real :: qsurf
    … …  
    1415      real :: zsurf
    1516      real :: albedo
     17      real :: snowmass
    1618
    1719      real :: time
    … …  
    3133
    3234      common/com_par1d/                                                 &
    33      & nat_surf,tsurf,rugos,                                            &
     35     & nat_surf,tsurf,rugos,rugosh,                                     &
    3436     & xqsol,qsurf,psurf,zsurf,albedo,time,time_ini,xlat,xlon,airefi,   &
    3537     & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp,            &
    3638     & forcing_type,tend_u,tend_v,tend_w,tend_t,tend_q,tend_rayo,       &
    3739     & nudge_u,nudge_v,nudge_w,nudge_t,nudge_q,                         &
    38      & iflag_nudge,                                                     &
     40     & iflag_nudge,snowmass,                                            &
    3941     & restart,ok_old_disvert
    4042

  • LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90

    r2635 r2672  
    2121       zgam, zmax0, zmea, zpic, zsig, &
    2222       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl
     23 
    2324   USE dimphy
    2425   USE surface_data, only : type_ocean,ok_veget
    … …  
    131132        logical :: forcing_amma    = .false.
    132133        logical :: forcing_dice    = .false.
     134        logical :: forcing_gabls4  = .false.
     135
    133136        logical :: forcing_GCM2SCM = .false.
    134137        logical :: forcing_GCSSold = .false.
    … …  
    174177      real :: pzero=1.e5
    175178      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),poub
     179      real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1)
    177180
    178181!---------------------------------------------------------------------
    … …  
    322325!             Different stages: soil model alone, atm. model alone
    323326!             then both models coupled
     327!forcing_type = 8 ==> forcing_gabls4 = .true.
     328!             initial profiles and large scale forcings in gabls4_driver.nc
    324329!forcing_type >= 100 ==> forcing_case = .true.
    325330!             initial profiles and large scale forcings in cas.nc
    … …  
    363368      elseif (forcing_type .eq.7) THEN
    364369       forcing_dice = .true.
     370      elseif (forcing_type .eq.8) THEN
     371       forcing_gabls4 = .true.
    365372      elseif (forcing_type .eq.101) THEN ! Cindynamo starts 1-10-2011 0h
    366373       forcing_case = .true.
    … …  
    457464      IF(forcing_type .EQ. 6) fnday=64800./86400.
    458465!     IF(forcing_type .EQ. 6) fnday=50400./86400.
     466 IF(forcing_type .EQ. 8 ) fnday=129600./86400.
    459467      annee_ref = anneeref
    460468      mois = 1
    … …  
    487495     & (year_ini_dice,mth_ini_dice,day_ini_dice,heure_ini_dice             &
    488496     & ,day_ju_ini_dice)
     497 ELSEIF (forcing_type .eq.8 ) THEN
     498! Convert the initial date of GABLS4 to Julian day
     499      call ymds2ju                                                         &
     500     & (year_ini_gabls4,mth_ini_gabls4,day_ini_gabls4,heure_ini_gabls4     &
     501     & ,day_ju_ini_gabls4)
    489502      ELSEIF (forcing_type .gt.100) THEN
    490503! Convert the initial date to Julian day
    … …  
    699712
    700713        fder=0.
    701         snsrf(1,:)=0.        ! couverture de neige des sous surface
     714        snsrf(1,:)=snowmass ! masse de neige des sous surface
    702715        qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
    703716        fevap=0.
    704717        z0m(1,:)=rugos     ! couverture de neige des sous surface
    705         z0h(1,:)=rugos     ! couverture de neige des sous surface
     718        z0h(1,:)=rugosh    ! couverture de neige des sous surface
    706719        agesno  = xagesno
    707720        tsoil(:,:,:)=tsurf
    … …  
    726739        print*,'avant phyredem'
    727740        pctsrf(1,:)=0.
    728         if (nat_surf.eq.0.) then
     741          if (nat_surf.eq.0.) then
    729742          pctsrf(1,is_oce)=1.
    730743          pctsrf(1,is_ter)=0.
    731         else
     744          pctsrf(1,is_lic)=0.
     745          pctsrf(1,is_sic)=0.
     746        else if (nat_surf .eq. 1) then
    732747          pctsrf(1,is_oce)=0.
    733748          pctsrf(1,is_ter)=1.
    734         end if
     749          pctsrf(1,is_lic)=0.
     750          pctsrf(1,is_sic)=0.
     751        else if (nat_surf .eq. 2) then
     752          pctsrf(1,is_oce)=0.
     753          pctsrf(1,is_ter)=0.
     754          pctsrf(1,is_lic)=1.
     755          pctsrf(1,is_sic)=0.
     756        else if (nat_surf .eq. 3) then
     757          pctsrf(1,is_oce)=0.
     758          pctsrf(1,is_ter)=0.
     759          pctsrf(1,is_lic)=0.
     760          pctsrf(1,is_sic)=1.
     761
     762     end if
     763
    735764
    736765        print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
Note: See TracChangeset for help on using the changeset viewer.