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

Introduction du cas GABLS4 par Etienne Vignon

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.