Ignore:
Timestamp:
Nov 21, 2019, 4:43:45 PM (4 years ago)
Author:
lguez
Message:

Merge revisions 3427:3600 of trunk into branch Ocean_skin

Location:
LMDZ6/branches/Ocean_skin
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/Ocean_skin

  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/mod_1D_cases_read2.F90

    r2764 r3605  
    315315END SUBROUTINE read2_1D_cas
    316316
     317!**********************************************************************************************
     318SUBROUTINE read_SCM_cas
     319      implicit none
     320
     321#include "netcdf.inc"
     322#include "date_cas.h"
     323
     324      INTEGER nid,rid,ierr
     325      INTEGER ii,jj,timeid
     326      REAL, ALLOCATABLE :: time_val(:)
     327
     328      print*,'ON EST VRAIMENT LA'
     329      fich_cas='cas.nc'
     330      print*,'fich_cas ',fich_cas
     331      ierr = NF_OPEN(fich_cas,NF_NOWRITE,nid)
     332      print*,'fich_cas,NF_NOWRITE,nid ',fich_cas,NF_NOWRITE,nid
     333      if (ierr.NE.NF_NOERR) then
     334         write(*,*) 'ERROR: GROS Pb opening forcings nc file '
     335         write(*,*) NF_STRERROR(ierr)
     336         stop ""
     337      endif
     338!.......................................................................
     339      ierr=NF_INQ_DIMID(nid,'lat',rid)
     340      IF (ierr.NE.NF_NOERR) THEN
     341         print*, 'Oh probleme lecture dimension lat'
     342      ENDIF
     343      ierr=NF_INQ_DIMLEN(nid,rid,ii)
     344      print*,'OK1 read2: nid,rid,lat',nid,rid,ii
     345!.......................................................................
     346      ierr=NF_INQ_DIMID(nid,'lon',rid)
     347      IF (ierr.NE.NF_NOERR) THEN
     348         print*, 'Oh probleme lecture dimension lon'
     349      ENDIF
     350      ierr=NF_INQ_DIMLEN(nid,rid,jj)
     351      print*,'OK2 read2: nid,rid,lat',nid,rid,jj
     352!.......................................................................
     353      ierr=NF_INQ_DIMID(nid,'lev',rid)
     354      IF (ierr.NE.NF_NOERR) THEN
     355         print*, 'Oh probleme lecture dimension nlev'
     356      ENDIF
     357      ierr=NF_INQ_DIMLEN(nid,rid,nlev_cas)
     358      print*,'OK3 read2: nid,rid,nlev_cas',nid,rid,nlev_cas
     359      IF ( .NOT. ( nlev_cas > 10 .AND. nlev_cas < 1000 )) THEN
     360              print*,'Valeur de nlev_cas peu probable'
     361              STOP
     362      ENDIF
     363!.......................................................................
     364      ierr=NF_INQ_DIMID(nid,'time',rid)
     365      nt_cas=0
     366      IF (ierr.NE.NF_NOERR) THEN
     367        stop 'Oh probleme lecture dimension time'
     368      ENDIF
     369      ierr=NF_INQ_DIMLEN(nid,rid,nt_cas)
     370      print*,'OK4 read2: nid,rid,nt_cas',nid,rid,nt_cas
     371! Lecture de l'axe des temps
     372      print*,'LECTURE DU TEMPS'
     373      ierr=NF_INQ_VARID(nid,'time',timeid)
     374         if(ierr/=NF_NOERR) then
     375           print *,'Variable time manquante dans cas.nc:'
     376           ierr=NF_NOERR
     377         else
     378                 allocate(time_val(nt_cas))
     379#ifdef NC_DOUBLE
     380         ierr = NF_GET_VAR_DOUBLE(nid,timeid,time_val)
     381#else
     382           ierr = NF_GET_VAR_REAL(nid,timeid,time_val)
     383#endif
     384           if(ierr/=NF_NOERR) then
     385              print *,'Pb a la lecture de time cas.nc: '
     386           endif
     387   endif
     388   IF (nt_cas>1) THEN
     389           pdt_cas=time_val(2)-time_val(1)
     390   ELSE
     391           pdt_cas=0.
     392   ENDIF
     393
     394
     395!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     396!profils moyens:
     397        allocate(plev_cas(nlev_cas,nt_cas),plevh_cas(nlev_cas+1))       
     398        allocate(z_cas(nlev_cas,nt_cas),zh_cas(nlev_cas+1))
     399        allocate(ap_cas(nlev_cas+1),bp_cas(nt_cas+1))
     400        allocate(t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas),qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas), &
     401             qi_cas(nlev_cas,nt_cas),rh_cas(nlev_cas,nt_cas))
     402        allocate(th_cas(nlev_cas,nt_cas),thl_cas(nlev_cas,nt_cas),thv_cas(nlev_cas,nt_cas),rv_cas(nlev_cas,nt_cas))
     403        allocate(u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas),vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas))
     404
     405!forcing
     406        allocate(ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas),dt_cas(nlev_cas,nt_cas),dtrad_cas(nlev_cas,nt_cas))
     407        allocate(hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas),dq_cas(nlev_cas,nt_cas))
     408        allocate(hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas),dth_cas(nlev_cas,nt_cas))
     409        allocate(hr_cas(nlev_cas,nt_cas),vr_cas(nlev_cas,nt_cas),dr_cas(nlev_cas,nt_cas))
     410        allocate(hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas),du_cas(nlev_cas,nt_cas))
     411        allocate(hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas),dv_cas(nlev_cas,nt_cas))
     412        allocate(ug_cas(nlev_cas,nt_cas))
     413        allocate(vg_cas(nlev_cas,nt_cas))
     414        allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas),ps_cas(nt_cas),ustar_cas(nt_cas),tke_cas(nt_cas))
     415        allocate(uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas),q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas))
     416
     417
     418
     419!champs interpoles
     420        allocate(plev_prof_cas(nlev_cas))
     421        allocate(t_prof_cas(nlev_cas))
     422        allocate(theta_prof_cas(nlev_cas))
     423        allocate(thl_prof_cas(nlev_cas))
     424        allocate(thv_prof_cas(nlev_cas))
     425        allocate(q_prof_cas(nlev_cas))
     426        allocate(qv_prof_cas(nlev_cas))
     427        allocate(ql_prof_cas(nlev_cas))
     428        allocate(qi_prof_cas(nlev_cas))
     429        allocate(rh_prof_cas(nlev_cas))
     430        allocate(rv_prof_cas(nlev_cas))
     431        allocate(u_prof_cas(nlev_cas))
     432        allocate(v_prof_cas(nlev_cas))
     433        allocate(vitw_prof_cas(nlev_cas))
     434        allocate(omega_prof_cas(nlev_cas))
     435        allocate(ug_prof_cas(nlev_cas))
     436        allocate(vg_prof_cas(nlev_cas))
     437        allocate(ht_prof_cas(nlev_cas))
     438        allocate(hth_prof_cas(nlev_cas))
     439        allocate(hq_prof_cas(nlev_cas))
     440        allocate(hu_prof_cas(nlev_cas))
     441        allocate(hv_prof_cas(nlev_cas))
     442        allocate(vt_prof_cas(nlev_cas))
     443        allocate(vth_prof_cas(nlev_cas))
     444        allocate(vq_prof_cas(nlev_cas))
     445        allocate(vu_prof_cas(nlev_cas))
     446        allocate(vv_prof_cas(nlev_cas))
     447        allocate(dt_prof_cas(nlev_cas))
     448        allocate(dth_prof_cas(nlev_cas))
     449        allocate(dtrad_prof_cas(nlev_cas))
     450        allocate(dq_prof_cas(nlev_cas))
     451        allocate(du_prof_cas(nlev_cas))
     452        allocate(dv_prof_cas(nlev_cas))
     453        allocate(uw_prof_cas(nlev_cas))
     454        allocate(vw_prof_cas(nlev_cas))
     455        allocate(q1_prof_cas(nlev_cas))
     456        allocate(q2_prof_cas(nlev_cas))
     457
     458        print*,'Allocations OK'
     459        call read_SCM (nid,nlev_cas,nt_cas,                                                                     &
     460     &     ap_cas,bp_cas,z_cas,plev_cas,zh_cas,plevh_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,                    &
     461     &     ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,ug_cas,vg_cas,du_cas,hu_cas,vu_cas,        &
     462     &     dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,               &
     463     &     dr_cas,hr_cas,vr_cas,dtrad_cas,sens_cas,lat_cas,ts_cas,ps_cas,ustar_cas,tke_cas,                      &
     464     &     uw_cas,vw_cas,q1_cas,q2_cas,orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough, &
     465     &     o3_cas,rugos_cas,clay_cas,sand_cas)
     466        print*,'Read2 cas OK'
     467        do ii=1,nlev_cas
     468        print*,'apres read2_cas, plev_cas=',ii,plev_cas(ii,1)
     469        enddo
     470
     471
     472END SUBROUTINE read_SCM_cas
    317473
    318474
     
    685841!-----------------------------------------------------------------------
    686842
     843
    687844         return
    688845         end subroutine read2_cas
     846
     847!======================================================================
     848      subroutine read_SCM(nid,nlevel,ntime,                                       &
     849     &     ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,ug,vg,&
     850     &     du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq,                                    &
     851     &     dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,ps,ustar,tke,uw,vw,q1,q2,       &
     852     &     orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,          &
     853     &     heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas)
     854
     855!program reading forcing of the case study
     856      implicit none
     857#include "netcdf.inc"
     858
     859      integer ntime,nlevel,k,t
     860
     861      real ap(nlevel+1),bp(nlevel+1)
     862      real zz(nlevel,ntime),zzh(nlevel+1)
     863      real pp(nlevel,ntime),pph(nlevel+1)
     864!profils initiaux
     865      real temp0(nlevel),qv0(nlevel),ql0(nlevel),qi0(nlevel),u0(nlevel),v0(nlevel),tke0(nlevel)
     866      real pp0(nlevel)   
     867      real temp(nlevel,ntime),qv(nlevel,ntime),ql(nlevel,ntime),qi(nlevel,ntime),rh(nlevel,ntime)
     868      real theta(nlevel,ntime),thv(nlevel,ntime),thl(nlevel,ntime),rv(nlevel,ntime)
     869      real u(nlevel,ntime),v(nlevel,ntime),tke(nlevel,ntime)
     870      real ug(nlevel,ntime),vg(nlevel,ntime)
     871      real vitw(nlevel,ntime),omega(nlevel,ntime)
     872      real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
     873      real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
     874      real dt(nlevel,ntime),ht(nlevel,ntime),vt(nlevel,ntime)
     875      real dtrad(nlevel,ntime)
     876      real dq(nlevel,ntime),hq(nlevel,ntime),vq(nlevel,ntime)
     877      real dth(nlevel,ntime),hth(nlevel,ntime),vth(nlevel,ntime),hthl(nlevel,ntime)
     878      real dr(nlevel,ntime),hr(nlevel,ntime),vr(nlevel,ntime)
     879      real flat(ntime),sens(ntime),ustar(ntime)
     880      real uw(nlevel,ntime),vw(nlevel,ntime),q1(nlevel,ntime),q2(nlevel,ntime)
     881      real ts(ntime),ps(ntime)
     882      real orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas
     883      real apbp(nlevel+1),resul(nlevel,ntime),resul1(nlevel),resul2(ntime),resul3
     884
     885
     886      integer nid, ierr,ierr1,ierr2,rid,i
     887      integer nbvar3d
     888      parameter(nbvar3d=70)
     889      integer var3didin(nbvar3d),missing_var(nbvar3d)
     890      character*13 name_var(1:nbvar3d)
     891      data name_var/'coor_par_a','coor_par_b','height_h','pressure_h',&
     892     &'temp','qv','ql','qi','u','v','tke','pressure',&
     893     &'w','omega','ug','vg','uadv','uadvh','uadvv','vadv','vadvh','vadvv','tadv','tadvh','tadvv',&
     894     &'qvadv','qvadvh','qvadvv','thadv','thadvh','thadvv','thladvh','radv','radvh','radvv','radcool','q1','q2','ustress','vstress', &
     895     'rh',&
     896     &'height_f','pressure_forc','tempt','theta','thv','thl','qvt','qlt','qit','rv','ut','vt','tket',&
     897     &'sfc_sens_flx','sfc_lat_flx','ts','ps','ustar',&
     898     &'orog','albedo','emiss','t_skin','q_skin','mom_rough','heat_rough','o3','rugos','clay','sand'/
     899      do i=1,nbvar3d
     900        missing_var(i)=0.
     901      enddo
     902
     903!-----------------------------------------------------------------------
     904
     905     print*,'ON EST LA'
     906       do i=1,nbvar3d
     907         ierr=NF_INQ_VARID(nid,name_var(i),var3didin(i))
     908         if(ierr/=NF_NOERR) then
     909           print *,'Variable manquante dans cas.nc:',i,name_var(i)
     910           ierr=NF_NOERR
     911           missing_var(i)=1
     912         else
     913!-----------------------------------------------------------------------
     914           if(i.LE.4) then     ! Lecture des coord pression en (nlevelp1,lat,lon)
     915#ifdef NC_DOUBLE
     916           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),apbp)
     917#else
     918           ierr = NF_GET_VAR_REAL(nid,var3didin(i),apbp)
     919#endif
     920           print *,'read2_cas(apbp), on a lu ',i,name_var(i)
     921           if(ierr/=NF_NOERR) then
     922              print *,'Pb a la lecture de cas.nc: ',name_var(i)
     923              stop "getvarup"
     924           endif
     925!-----------------------------------------------------------------------
     926           else if(i.gt.4.and.i.LE.12) then   ! Lecture des variables en (time,nlevel,lat,lon)
     927#ifdef NC_DOUBLE
     928           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul1)
     929#else
     930           ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul1)
     931#endif
     932           print *,'read2_cas(resul1), on a lu ',i,name_var(i)
     933           if(ierr/=NF_NOERR) then
     934              print *,'Pb a la lecture de cas.nc: ',name_var(i)
     935              stop "getvarup"
     936           endif
     937         print*,'Lecture de la variable #i ',i,name_var(i),minval(resul1),maxval(resul1)
     938!-----------------------------------------------------------------------
     939           else if(i.gt.12.and.i.LE.54) then   ! Lecture des variables en (time,nlevel,lat,lon)
     940#ifdef NC_DOUBLE
     941           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul)
     942#else
     943           ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul)
     944#endif
     945           print *,'read2_cas(resul), on a lu ',i,name_var(i)
     946           if(ierr/=NF_NOERR) then
     947              print *,'Pb a la lecture de cas.nc: ',name_var(i)
     948              stop "getvarup"
     949           endif
     950         print*,'Lecture de la variable #i ',i,name_var(i),minval(resul),maxval(resul)
     951!-----------------------------------------------------------------------
     952           else if (i.gt.54.and.i.LE.65) then   ! Lecture des variables en (time,lat,lon)
     953#ifdef NC_DOUBLE
     954           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul2)
     955#else
     956           ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul2)
     957#endif
     958           print *,'read2_cas(resul2), on a lu ',i,name_var(i)
     959           if(ierr/=NF_NOERR) then
     960              print *,'Pb a la lecture de cas.nc: ',name_var(i)
     961              stop "getvarup"
     962           endif
     963         print*,'Lecture de la variable #i  ',i,name_var(i),minval(resul2),maxval(resul2)
     964!-----------------------------------------------------------------------
     965           else     ! Lecture des constantes (lat,lon)
     966#ifdef NC_DOUBLE
     967           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul3)
     968#else
     969           ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul3)
     970#endif
     971           print *,'read2_cas(resul3), on a lu ',i,name_var(i)
     972           if(ierr/=NF_NOERR) then
     973              print *,'Pb a la lecture de cas.nc: ',name_var(i)
     974              stop "getvarup"
     975           endif
     976         print*,'Lecture de la variable #i ',i,name_var(i),resul3
     977           endif
     978         endif
     979!-----------------------------------------------------------------------
     980         select case(i)
     981         !case(1) ; ap=apbp       ! donnees indexees en nlevel+1
     982         ! case(2) ; bp=apbp
     983           case(3) ; zzh=apbp
     984           case(4) ; pph=apbp
     985           case(5) ; temp0=resul1    ! donnees initiales
     986           case(6) ; qv0=resul1
     987           case(7) ; ql0=resul1
     988           case(8) ; qi0=resul1
     989           case(9) ; u0=resul1
     990           case(10) ; v0=resul1
     991           case(11) ; tke0=resul1
     992           case(12) ; pp0=resul1
     993           case(13) ; vitw=resul    ! donnees indexees en nlevel,time
     994           case(14) ; omega=resul
     995           case(15) ; ug=resul
     996           case(16) ; vg=resul
     997           case(17) ; du=resul
     998           case(18) ; hu=resul
     999           case(19) ; vu=resul
     1000           case(20) ; dv=resul
     1001           case(21) ; hv=resul
     1002           case(22) ; vv=resul
     1003           case(23) ; dt=resul
     1004           case(24) ; ht=resul
     1005           case(25) ; vt=resul
     1006           case(26) ; dq=resul
     1007           case(27) ; hq=resul
     1008           case(28) ; vq=resul
     1009           case(29) ; dth=resul
     1010           case(30) ; hth=resul
     1011           case(31) ; vth=resul
     1012           case(32) ; hthl=resul
     1013           case(33) ; dr=resul
     1014           case(34) ; hr=resul
     1015           case(35) ; vr=resul
     1016           case(36) ; dtrad=resul
     1017           case(37) ; q1=resul
     1018           case(38) ; q2=resul
     1019           case(39) ; uw=resul
     1020           case(40) ; vw=resul
     1021           case(41) ; rh=resul
     1022           case(42) ; zz=resul      ! donnees en time,nlevel pour profil initial
     1023           case(43) ; pp=resul
     1024           case(44) ; temp=resul
     1025           case(45) ; theta=resul
     1026           case(46) ; thv=resul
     1027           case(47) ; thl=resul
     1028           case(48) ; qv=resul
     1029           case(49) ; ql=resul
     1030           case(50) ; qi=resul
     1031           case(51) ; rv=resul
     1032           case(52) ; u=resul
     1033           case(53) ; v=resul
     1034           case(54) ; tke=resul
     1035           case(55) ; sens=resul2   ! donnees indexees en time
     1036           case(56) ; flat=resul2
     1037           case(57) ; ts=resul2
     1038           case(58) ; ps=resul2
     1039           case(59) ; ustar=resul2
     1040           case(60) ; orog_cas=resul3      ! constantes
     1041           case(61) ; albedo_cas=resul3
     1042           case(62) ; emiss_cas=resul3
     1043           case(63) ; t_skin_cas=resul3
     1044           case(64) ; q_skin_cas=resul3
     1045           case(65) ; mom_rough=resul3
     1046           case(66) ; heat_rough=resul3
     1047           case(67) ; o3_cas=resul3       
     1048           case(68) ; rugos_cas=resul3
     1049           case(69) ; clay_cas=resul3
     1050           case(70) ; sand_cas=resul3
     1051         end select
     1052         resul=0.
     1053         resul1=0.
     1054         resul2=0.
     1055         resul3=0.
     1056       enddo
     1057         print*,'Lecture de la variable APRES ,sens ',minval(sens),maxval(sens)
     1058         print*,'Lecture de la variable APRES ,flat ',minval(flat),maxval(flat)
     1059
     1060!CR:ATTENTION EN ATTENTE DE REGLER LA QUESTION DU PAS DE TEMPS INITIAL
     1061       do t=1,ntime
     1062          do k=1,nlevel
     1063             temp(k,t)=temp0(k)
     1064             qv(k,t)=qv0(k)
     1065             ql(k,t)=ql0(k)
     1066             qi(k,t)=qi0(k)
     1067             u(k,t)=u0(k)
     1068             v(k,t)=v0(k)
     1069             tke(k,t)=tke0(k)
     1070          enddo
     1071       enddo
     1072!-----------------------------------------------------------------------
     1073
     1074         return
     1075         end subroutine read_SCM
     1076!======================================================================
     1077
    6891078!======================================================================
    6901079        SUBROUTINE interp_case_time2(day,day1,annee_ref                &
Note: See TracChangeset for help on using the changeset viewer.