Ignore:
Timestamp:
Jan 11, 2021, 11:24:08 PM (4 years ago)
Author:
lguez
Message:

Sync latest trunk changes to Ocean_skin

Location:
LMDZ6/branches/Ocean_skin
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/Ocean_skin

  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1DUTILS.h

    r3605 r3798  
    233233       CALL getin('ok_flux_surf',ok_flux_surf)
    234234
     235!Config  Key  = ok_forc_tsurf
     236!Config  Desc = forcage ou non par la Ts
     237!Config  Def  = false
     238!Config  Help = forcage ou non par la Ts
     239       ok_forc_tsurf=.false.
     240       CALL getin('ok_forc_tsurf',ok_forc_tsurf)
     241
    235242!Config  Key  = ok_prescr_ust
    236243!Config  Desc = ustar impose ou non
     
    239246       ok_prescr_ust = .false.
    240247       CALL getin('ok_prescr_ust',ok_prescr_ust)
     248
     249
     250!Config  Key  = ok_prescr_beta
     251!Config  Desc = betaevap impose ou non
     252!Config  Def  = false
     253!Config  Help = betaevap impose ou non
     254       ok_prescr_beta = .false.
     255       CALL getin('ok_prescr_beta',ok_prescr_beta)
    241256
    242257!Config  Key  = ok_old_disvert
     
    280295!Config  Desc = surface temperature
    281296!Config  Def  = 290.
    282 !Config  Help = not used if type_ts_forcing=1 in lmdz1d.F
     297!Config  Help = surface temperature
    283298       tsurf = 290.
    284299       CALL getin('tsurf',tsurf)
     
    297312       zsurf = 0.
    298313       CALL getin('zsurf',zsurf)
     314! EV pour accord avec format standard       
     315       CALL getin('zorog',zsurf)
     316
    299317
    300318!Config  Key  = rugos
     
    304322       rugos = 0.0001
    305323       CALL getin('rugos',rugos)
     324! FH/2020/04/08/confinement: Pour le nouveau format standard, la rugosite s'appelle z0
     325       CALL getin('z0',rugos)
    306326
    307327!Config  Key  = rugosh
     
    357377       qsolinp = 1.
    358378       CALL getin('qsolinp',qsolinp)
     379
     380
     381
     382!Config  Key  = betaevap
     383!Config  Desc = beta for actual evaporation when prescribed
     384!Config  Def  = 1.0
     385!Config  Help =
     386       betaevap = 1.
     387       CALL getin('betaevap',betaevap)     
    359388
    360389!Config  Key  = zpicinp
     
    518547       CALL getin('forc_ustar',forc_ustar)
    519548       IF (forc_ustar .EQ. 1) ok_prescr_ust=.true.
     549
    520550
    521551!Config  Key  = nudging_u
     
    12461276      END
    12471277
    1248 !======================================================================
    1249        SUBROUTINE read_tsurf1d(knon,sst_out)
    1250 
    1251 ! This subroutine specifies the surface temperature to be used in 1D simulations
    1252 
    1253       USE dimphy, ONLY : klon
    1254 
    1255       INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
    1256       REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
    1257 
    1258        INTEGER :: i
    1259 ! COMMON defined in lmdz1d.F:
    1260        real ts_cur
    1261        common /sst_forcing/ts_cur
    1262 
    1263        DO i = 1, knon
    1264         sst_out(i) = ts_cur
    1265        ENDDO
    1266 
    1267       END SUBROUTINE read_tsurf1d
    1268 
     1278!!======================================================================
     1279!       SUBROUTINE read_tsurf1d(knon,sst_out)
     1280!
     1281!! This subroutine specifies the surface temperature to be used in 1D simulations
     1282!
     1283!      USE dimphy, ONLY : klon
     1284!
     1285!      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
     1286!      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
     1287!
     1288!       INTEGER :: i
     1289!! COMMON defined in lmdz1d.F:
     1290!       real ts_cur
     1291!       common /sst_forcing/ts_cur
     1292
     1293!       DO i = 1, knon
     1294!        sst_out(i) = ts_cur
     1295!       ENDDO
     1296!
     1297!      END SUBROUTINE read_tsurf1d
     1298!
    12691299!===============================================================
    12701300      subroutine advect_vert(llm,w,dt,q,plev)
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_decl_cases.h

    r3605 r3798  
    3434        real w_mod(llm), t_mod(llm),q_mod(llm)
    3535        real u_mod(llm),v_mod(llm), ht_mod(llm),vt_mod(llm),ug_mod(llm),vg_mod(llm)
    36         real temp_nudg_mod(llm),qv_nudg_mod(llm),u_nudg_mod(llm),v_nudg_mod(llm)
     36            real temp_nudg_mod(llm),qv_nudg_mod(llm),u_nudg_mod(llm),v_nudg_mod(llm)
    3737        real hq_mod(llm),vq_mod(llm),qv_mod(llm),ql_mod(llm),qt_mod(llm)
    3838        real th_mod(llm)
    3939
    40         real ts_cur
    41         common /sst_forcing/ts_cur ! also in read_tsurf1d.F
     40! EV comment these lines
     41!        real ts_cur
     42!        common /sst_forcing/ts_cur ! also in read_tsurf1d.F
    4243!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4344! Declarations specifiques au cas RICO
     
    286287        real temp_nudg_mod_cas(llm),qv_nudg_mod_cas(llm),v_nudg_mod_cas(llm),u_nudg_mod_cas(llm)
    287288        real u_mod_cas(llm),v_mod_cas(llm)
    288         real omega_mod_cas(llm)
     289        real omega_mod_cas(llm),tke_mod_cas(llm+1)
    289290        real ht_mod_cas(llm),vt_mod_cas(llm),dt_mod_cas(llm),dtrad_mod_cas(llm)
    290291        real hth_mod_cas(llm),vth_mod_cas(llm),dth_mod_cas(llm)
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_interp_cases.h

    r3605 r3798  
    11
    2          print*,'FORCING CASE forcing_case2'
     2        print*,'FORCING CASE forcing_case2'
    33!       print*,                                                             &
    44!    & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
     
    1313     &       ,u_cas,v_cas,ug_cas,vg_cas                                                     &
    1414     &       ,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas                               &
    15      &       ,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas                                       &
     15     &       ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas                                       &
    1616     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
    1717     &       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas       &
    18      &       ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                                           &
     18     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                                           &
    1919!
    2020     &       ,ts_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas  &
     
    2222     &       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                                 &
    2323     &       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas           &
    24      &       ,vitw_prof_cas,omega_prof_cas                                                  &
     24     &       ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                     &
    2525     &       ,du_prof_cas,hu_prof_cas,vu_prof_cas                                           &
    2626     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
    2727     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
    2828     &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
    29      &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
    30 
    31              ts_cur = ts_prof_cas
     29     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
     30! EV tg instead of ts_cur
     31             tg = ts_prof_cas
    3232!            psurf=plev_prof_cas(1)
    3333             psurf=ps_prof_cas
    3434
    3535! vertical interpolation:
    36       CALL interp2_case_vertical_std(play,nlev_cas,plev_prof_cas                                              &
     36      CALL interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                                              &
    3737     &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                       &
    3838     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
    3939     &         ,ug_prof_cas,vg_prof_cas                                                                   &
    4040     &         ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                       &
    41      &         ,vitw_prof_cas,omega_prof_cas                                      &
     41     &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                                 &
    4242     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
    4343     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
     
    4747     &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas                                                 &
    4848     &         ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                           &
    49      &         ,w_mod_cas,omega_mod_cas                                                                   &
     49     &         ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                       &
    5050     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
    5151     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
     
    109109      do l = 1, llm
    110110! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
     111        print*, l, llm
     112        print*, play(l), temp(l)
    111113       omega(l) = -w_mod_cas(l)*play(l)*rg/(rd*temp(l))
    112114      enddo
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_read_forc_cases.h

    r3605 r3798  
    2727     &       ,u_cas,v_cas,ug_cas,vg_cas                                                     &
    2828     &       ,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas                               &
    29      &       ,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas                                       &
     29     &       ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas                                       &
    3030     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
    3131     &       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas       &
    32      &       ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                                           &
     32     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                                           &
    3333!
    3434     &       ,ts_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas  &
     
    3636     &       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                                 &
    3737     &       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas           &
    38      &       ,vitw_prof_cas,omega_prof_cas                                                  &
     38     &       ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                    &
    3939     &       ,du_prof_cas,hu_prof_cas,vu_prof_cas                                           &
    4040     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
    4141     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
    4242     &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
    43      &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
     43     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
    4444
    4545      do l = 1, nlev_cas
     
    4949! vertical interpolation using interpolation routine:
    5050!      write(*,*)'avant interp vert', t_prof
    51       CALL interp2_case_vertical_std(play,nlev_cas,plev_prof_cas                                              &
     51      CALL interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                                              &
    5252     &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                          &
    5353     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
     
    5555     &       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                         &
    5656
    57      &         ,vitw_prof_cas,omega_prof_cas                                                              &
     57     &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                                 &
    5858     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
    5959     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
     
    6363     &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas                                                 &
    6464     &         ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                           &
    65      &         ,w_mod_cas,omega_mod_cas                                                                   &
     65     &         ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                       &
    6666     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
    6767     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
     
    7070
    7171! initial and boundary conditions :
    72 !      tsurf = ts_prof_cas
     72!     tsurf = ts_prof_cas
    7373      psurf = ps_prof_cas
    74       ts_cur = ts_prof_cas
     74      !EV tg instead of ts_cur
     75      tg = ts_prof_cas
     76      print*, 'tg=', tg
     77
    7578      do l = 1, llm
    7679       temp(l) = t_mod_cas(l)
     
    9598       d_u_adv(l) = du_mod_cas(l)+hu_mod_cas(l)+vu_mod_cas(l)
    9699       d_v_adv(l) = dv_mod_cas(l)+hv_mod_cas(l)+vv_mod_cas(l)
     100      enddo
    97101
    98 !print*,'d_t_adv ',d_t_adv(1:20)*86400
     102! Etienne pour initialisation de TKE
    99103
    100       enddo     
     104       do l=1,llm+1
     105       pbl_tke(:,l,:)=tke_mod_cas(l)
     106       enddo     
    101107
    102108! Faut-il multiplier par -1 ? (MPL 20160713)
     
    108114       IF (ok_prescr_ust) THEN
    109115       ust=ustar_prof_cas
    110        print *,'ust=',ust
    111116       ENDIF
    112117
     118
    113119      endif !forcing_SCM
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/mod_1D_cases_read2.F90

    r3605 r3798  
    316316
    317317!**********************************************************************************************
    318 SUBROUTINE read_SCM_cas
     318SUBROUTINE old_read_SCM_cas
    319319      implicit none
    320320
     
    457457
    458458        print*,'Allocations OK'
    459         call read_SCM (nid,nlev_cas,nt_cas,                                                                     &
     459        call old_read_SCM (nid,nlev_cas,nt_cas,                                                                     &
    460460     &     ap_cas,bp_cas,z_cas,plev_cas,zh_cas,plevh_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,                    &
    461461     &     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,        &
     
    470470
    471471
    472 END SUBROUTINE read_SCM_cas
     472END SUBROUTINE old_read_SCM_cas
    473473
    474474
     
    846846
    847847!======================================================================
    848       subroutine read_SCM(nid,nlevel,ntime,                                       &
     848      subroutine old_read_SCM(nid,nlevel,ntime,                                       &
    849849     &     ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,ug,vg,&
    850850     &     du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq,                                    &
     
    10731073
    10741074         return
    1075          end subroutine read_SCM
     1075         end subroutine old_read_SCM
    10761076!======================================================================
    10771077
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90

    r3605 r3798  
    1818        real, allocatable::  t_cas(:,:),q_cas(:,:),qv_cas(:,:),ql_cas(:,:),qi_cas(:,:),rh_cas(:,:)
    1919        real, allocatable::  th_cas(:,:),thv_cas(:,:),thl_cas(:,:),rv_cas(:,:)
    20         real, allocatable::  u_cas(:,:),v_cas(:,:),vitw_cas(:,:),omega_cas(:,:)
     20        real, allocatable::  u_cas(:,:),v_cas(:,:),vitw_cas(:,:),omega_cas(:,:),tke_cas(:,:)
    2121
    2222!forcing
     
    3030        real, allocatable::  temp_nudg_cas(:,:),qv_nudg_cas(:,:),u_nudg_cas(:,:),v_nudg_cas(:,:)
    3131        real, allocatable::  lat_cas(:),sens_cas(:),ts_cas(:),ps_cas(:),ustar_cas(:)
    32         real, allocatable::  uw_cas(:,:),vw_cas(:,:),q1_cas(:,:),q2_cas(:,:),tke_cas(:)
     32        real, allocatable::  uw_cas(:,:),vw_cas(:,:),q1_cas(:,:),q2_cas(:,:),tkes_cas(:)
    3333
    3434!champs interpoles
     
    4848        real, allocatable::  vitw_prof_cas(:)
    4949        real, allocatable::  omega_prof_cas(:)
     50        real, allocatable::  tke_prof_cas(:)
    5051        real, allocatable::  ug_prof_cas(:)
    5152        real, allocatable::  vg_prof_cas(:)
     
    7374
    7475
    75         real lat_prof_cas,sens_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas,tke_prof_cas
     76        real lat_prof_cas,sens_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas,tkes_prof_cas
    7677        real o3_cas,orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough,rugos_cas,sand_cas,clay_cas
    7778     
     
    9293      REAL, ALLOCATABLE :: time_val(:)
    9394
    94       print*,'ON EST VRAIMENT LA'
     95      print*,'ON EST VRAIMENT DASN MOD_1D_CASES_READ_STD'
    9596      fich_cas='cas.nc'
    9697      print*,'fich_cas ',fich_cas
     
    123124      ierr=NF_INQ_DIMLEN(nid,rid,nlev_cas)
    124125      print*,'OK3 read2: nid,rid,nlev_cas',nid,rid,nlev_cas
    125       IF ( .NOT. ( nlev_cas > 10 .AND. nlev_cas < 1000 )) THEN
     126      IF ( .NOT. ( nlev_cas > 10 .AND. nlev_cas < 200000 )) THEN
    126127              print*,'Valeur de nlev_cas peu probable'
    127128              STOP
     
    168169        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))
    169170        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))
    170 
     171        allocate(tke_cas(nlev_cas,nt_cas))
    171172!forcing
    172173        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))
     
    179180        allocate(temp_nudg_cas(nlev_cas,nt_cas),qv_nudg_cas(nlev_cas,nt_cas))
    180181        allocate(u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas))
    181         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))
     182        allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas),ps_cas(nt_cas),ustar_cas(nt_cas),tkes_cas(nt_cas))
    182183        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))
    183184
     
    200201        allocate(vitw_prof_cas(nlev_cas))
    201202        allocate(omega_prof_cas(nlev_cas))
     203        allocate(tke_prof_cas(nlev_cas))
    202204        allocate(ug_prof_cas(nlev_cas))
    203205        allocate(vg_prof_cas(nlev_cas))
     
    228230        CALL read_SCM (nid,nlev_cas,nt_cas,                                                                     &
    229231     &     ap_cas,bp_cas,z_cas,plev_cas,zh_cas,plevh_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,                   &
    230      &     ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,ug_cas,vg_cas,                            &
     232     &     ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,tke_cas,ug_cas,vg_cas,                            &
    231233     &     temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas,                                                     &
    232234     &     du_cas,hu_cas,vu_cas,                                                                                &
    233235     &     dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,               &
    234      &     dr_cas,hr_cas,vr_cas,dtrad_cas,sens_cas,lat_cas,ts_cas,ps_cas,ustar_cas,tke_cas,                      &
     236     &     dr_cas,hr_cas,vr_cas,dtrad_cas,sens_cas,lat_cas,ts_cas,ps_cas,ustar_cas,tkes_cas,                      &
    235237     &     uw_cas,vw_cas,q1_cas,q2_cas,orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough, &
    236238     &     o3_cas,rugos_cas,clay_cas,sand_cas)
     
    254256        deallocate(t_cas,q_cas,qv_cas,ql_cas,qi_cas,rh_cas)
    255257        deallocate(th_cas,thl_cas,thv_cas,rv_cas)
    256         deallocate(u_cas,v_cas,vitw_cas,omega_cas)
     258        deallocate(u_cas,v_cas,vitw_cas,omega_cas,tke_cas)
    257259       
    258260!forcing
     
    265267        deallocate(ug_cas)
    266268        deallocate(vg_cas)
    267         deallocate(lat_cas,sens_cas,ts_cas,ps_cas,ustar_cas,tke_cas,uw_cas,vw_cas,q1_cas,q2_cas)
     269        deallocate(lat_cas,sens_cas,ts_cas,ps_cas,ustar_cas,tkes_cas,uw_cas,vw_cas,q1_cas,q2_cas)
    268270
    269271!champs interpoles
     
    283285        deallocate(vitw_prof_cas)
    284286        deallocate(omega_prof_cas)
     287        deallocate(tke_prof_cas)
    285288        deallocate(ug_prof_cas)
    286289        deallocate(vg_prof_cas)
     
    312315!=====================================================================
    313316      SUBROUTINE read_SCM(nid,nlevel,ntime,                                       &
    314      &     ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,ug,vg,&
     317     &     ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,tke,ug,vg,&
    315318     &     temp_nudg,qv_nudg,u_nudg,v_nudg,                                        &
    316319     &     du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq,                                    &
    317      &     dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,ps,ustar,tke,uw,vw,q1,q2,       &
     320     &     dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,ps,ustar,tkes,uw,vw,q1,q2,       &
    318321     &     orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,          &
    319322     &     heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas)
     
    334337      real temp(nlevel,ntime),qv(nlevel,ntime),ql(nlevel,ntime),qi(nlevel,ntime),rh(nlevel,ntime)
    335338      real theta(nlevel,ntime),thv(nlevel,ntime),thl(nlevel,ntime),rv(nlevel,ntime)
    336       real u(nlevel,ntime),v(nlevel,ntime),tke(nlevel,ntime)
     339      real u(nlevel,ntime),v(nlevel,ntime),tkes(ntime)
    337340      real temp_nudg(nlevel,ntime),qv_nudg(nlevel,ntime),u_nudg(nlevel,ntime),v_nudg(nlevel,ntime)
    338341      real ug(nlevel,ntime),vg(nlevel,ntime)
    339       real vitw(nlevel,ntime),omega(nlevel,ntime)
     342      real vitw(nlevel,ntime),omega(nlevel,ntime),tke(nlevel,ntime)
    340343      real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
    341344      real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
     
    365368     &'temp','qv','ql','qi','u','v','tke','pressure',& ! #5-#12
    366369     ! coordonnees pression + temps #42
    367      &'w','omega','ug','vg','uadv','uadvh','uadvv','vadv','vadvh','vadvv','tadv','tadvh','tadvv',& !  #13 - #25
    368      &'qvadv','qvadvh','qvadvv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #33
    369      & 'radv','radvh','radvv','radcool','q1','q2','ustress','vstress',                           & ! #34 - #41
    370      & 'rh','temp_nudg','qv_nudg','u_nudg','v_nudg',&
    371      &'height_f','pressure_forc','tempt','theta','thv','thl','qvt','qlt','qit','rv','ut','vt','tket',&
     370     &'w','omega','ug','vg','uadv','uadvh','uadvv','vadv','vadvh','vadvv','temp_adv','tadvh','tadvv',& !  #13 - #25
     371     &'qv_adv','qvadvh','qvadvv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #32
     372     & 'radv','radvh','radvv','radcool','q1','q2','ustress','vstress',                           & ! #33 - #40
     373     & 'rh','temp_nudging','qv_nudging','u_nudging','v_nudging',                                       & ! #41-45
     374     &'height_f','pressure_forc','tempt','theta','thv','thl','qvt','qlt','qit','rv','ut','vt',   & ! #46-58
    372375     ! coordonnees temps #12
    373      &'sfc_sens_flx','sfc_lat_flx','ts','ps','ustar',&
     376     &'tkes','sfc_sens_flx','sfc_lat_flx','ts','ps','ustar',&
    374377     &'orog','albedo','emiss','t_skin','q_skin','mom_rough','heat_rough',&
    375378     ! scalaires #4
    376379     &'o3','rugos','clay','sand'/
    377380
    378       do i=1,nbvar3d
    379         missing_var(i)=0.
    380       enddo
    381 
     381!-----------------------------------------------------------------------
     382! Checking availability of variable #i in the cas.nc file
     383!     missing_var=1 if the variable is missing
    382384!-----------------------------------------------------------------------
    383385
    384386       do i=1,nbvar3d
     387         missing_var(i)=0.
    385388         ierr=NF_INQ_VARID(nid,name_var(i),var3didin(i))
    386389         if(ierr/=NF_NOERR) then
     
    391394
    392395!-----------------------------------------------------------------------
    393 ! Activation de quelques cles en fonction des variables disponibles
     396! Activating keys depending on the presence of specific variables in cas.nc
    394397!-----------------------------------------------------------------------
    395398if ( 1 == 1 ) THEN
    396             if ( name_var(i) == 'temp_nudg' .and. nint(nudging_t)==0) stop 'Nudging inconsistency temp'
    397             if ( name_var(i) == 'qv_nudg' .and. nint(nudging_qv)==0) stop 'Nudging inconsistency qv'
    398             if ( name_var(i) == 'u_nudg' .and. nint(nudging_u)==0) stop 'Nudging inconsistency u'
    399             if ( name_var(i) == 'v_nudg' .and. nint(nudging_u)==0) stop 'Nudging inconsistency v'
     399! A MODIFIER: il faudrait dire nudging_temp mais faut le declarer dans compar1d.h etc...       
     400!           if ( name_var(i) == 'temp_nudging' .and. nint(nudging_t)==0) stop 'Nudging inconsistency temp'
     401            if ( name_var(i) == 'qv_nudging' .and. nint(nudging_qv)==0) stop 'Nudging inconsistency qv'
     402            if ( name_var(i) == 'u_nudging' .and. nint(nudging_u)==0) stop 'Nudging inconsistency u'
     403            if ( name_var(i) == 'v_nudging' .and. nint(nudging_u)==0) stop 'Nudging inconsistency v'
    400404    ELSE
    401405             print*,'GUIDAGE : CONSISTENCY CHECK DEACTIVATED FOR TESTS of SANDU/REF'
     
    403407
    404408!-----------------------------------------------------------------------
    405            if(i.LE.4) then     ! Lecture des coord pression en (nlevelp1,lat,lon)
     409! Reading variables 1D (N+1) vertical variables (nlevelp1,lat,lon)
     410!-----------------------------------------------------------------------
     411           if(i.LE.4) then
    406412#ifdef NC_DOUBLE
    407413           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),apbp)
     
    414420              stop "getvarup"
    415421           endif
    416 !-----------------------------------------------------------------------
    417            else if(i.gt.4.and.i.LE.12) then   ! Lecture des variables en (time,nlevel,lat,lon)
     422
     423!-----------------------------------------------------------------------
     424!  Reading 1D (N) vertical varialbes    (nlevel,lat,lon)   
     425!-----------------------------------------------------------------------
     426           else if(i.gt.4.and.i.LE.12) then 
    418427#ifdef NC_DOUBLE
    419428           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul1)
     
    427436           endif
    428437         print*,'Lecture de la variable #i ',i,name_var(i),minval(resul1),maxval(resul1)
    429 !-----------------------------------------------------------------------
    430            else if(i.gt.12.and.i.LE.54) then   ! Lecture des variables en (time,nlevel,lat,lon)
     438
     439!-----------------------------------------------------------------------
     440!  Reading 2D tim-vertical variables  (time,nlevel,lat,lon)
     441!  TBD : seems to be the same as above.
     442!-----------------------------------------------------------------------
     443           else if(i.gt.12.and.i.LE.57) then
    431444#ifdef NC_DOUBLE
    432445           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul)
     
    439452              stop "getvarup"
    440453           endif
    441 
    442454         print*,'Lecture de la variable #i ',i,name_var(i),minval(resul),maxval(resul)
    443 !-----------------------------------------------------------------------
    444            else if (i.gt.54.and.i.LE.65) then   ! Lecture des variables en (time,lat,lon)
     455
     456!-----------------------------------------------------------------------
     457!  Reading 1D time variables (time,lat,lon)
     458!-----------------------------------------------------------------------
     459           else if (i.gt.57.and.i.LE.63) then
    445460#ifdef NC_DOUBLE
    446461           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul2)
     
    454469           endif
    455470         print*,'Lecture de la variable #i  ',i,name_var(i),minval(resul2),maxval(resul2)
    456 !-----------------------------------------------------------------------
    457            else     ! Lecture des constantes (lat,lon)
     471
     472!-----------------------------------------------------------------------
     473! Reading scalar variables (lat,lon)
     474!-----------------------------------------------------------------------
     475           else
    458476#ifdef NC_DOUBLE
    459477           ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul3)
     
    469487           endif
    470488         endif
     489
     490!-----------------------------------------------------------------------
     491! Attributing variables
    471492!-----------------------------------------------------------------------
    472493         select case(i)
     
    528549           case(56) ; u=resul
    529550           case(57) ; v=resul
    530            case(58) ; tke=resul
    531            case(59) ; sens=resul2   ! donnees indexees en time
     551           case(58) ; tkes=resul2   ! donnees indexees en time
     552           case(59) ; sens=resul2
    532553           case(60) ; flat=resul2
    533554           case(61) ; ts=resul2
     
    577598
    578599!**********************************************************************************************
     600
     601!**********************************************************************************************
    579602        SUBROUTINE interp_case_time_std(day,day1,annee_ref                           &
    580603!    &         ,year_cas,day_cas,nt_cas,pdt_forc,nlev_cas                         &
     
    583606     &         ,qv_cas,ql_cas,qi_cas,u_cas,v_cas                                  &
    584607     &         ,ug_cas,vg_cas,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas     &
    585      &         ,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas             &
     608     &         ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas             &
    586609     &         ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas               &
    587610     &         ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas                      &
    588611     &         ,lat_cas,sens_cas,ustar_cas                                        &
    589      &         ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                               &
     612     &         ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                               &
    590613!
    591614     &         ,ts_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas               &
     
    593616     &         ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                     &
    594617     &         ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas     &
    595      &         ,vitw_prof_cas,omega_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas  &
     618     &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas  &
    596619     &         ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas                   &
    597620     &         ,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas                &
    598621     &         ,hq_prof_cas,vq_prof_cas,dth_prof_cas,hth_prof_cas,vth_prof_cas    &
    599622     &         ,lat_prof_cas,sens_prof_cas                                        &
    600      &         ,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
    601          
    602 
    603         implicit none
     623     &         ,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
     624
     625
     626
     627
     628
     629
     630       implicit none
    604631
    605632!---------------------------------------------------------------------------------------
     
    621648        real ts_cas(nt_cas),ps_cas(nt_cas)
    622649        real plev_cas(nlev_cas,nt_cas)
    623         real t_cas(nlev_cas,nt_cas),theta_cas(nlev_cas,nt_cas),thv_cas(nlev_cas,nt_cas),thl_cas(nlev_cas,nt_cas)
     650        real t_cas(nlev_cas,nt_cas),theta_cas(nlev_cas,nt_cas)
     651        real thv_cas(nlev_cas,nt_cas), thl_cas(nlev_cas,nt_cas)
    624652        real qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas),qi_cas(nlev_cas,nt_cas)
    625653        real u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas)
     
    628656        real u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas)
    629657
    630         real vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas)
     658        real vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas),tke_cas(nlev_cas,nt_cas)
    631659        real du_cas(nlev_cas,nt_cas),hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas)
    632660        real dv_cas(nlev_cas,nt_cas),hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas)
     
    635663        real dtrad_cas(nlev_cas,nt_cas)
    636664        real dq_cas(nlev_cas,nt_cas),hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas)
    637         real lat_cas(nt_cas),sens_cas(nt_cas),tke_cas(nt_cas)
     665        real lat_cas(nt_cas),sens_cas(nt_cas),tkes_cas(nt_cas)
    638666        real ustar_cas(nt_cas),uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas)
    639667        real q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas)
     
    648676        real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
    649677
    650         real vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
     678        real vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas),tke_prof_cas(nlev_cas)
    651679        real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
    652680        real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
     
    655683        real dtrad_prof_cas(nlev_cas)
    656684        real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
    657         real lat_prof_cas,sens_prof_cas,tke_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas
     685        real lat_prof_cas,sens_prof_cas,tkes_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas
    658686        real uw_prof_cas(nlev_cas),vw_prof_cas(nlev_cas),q1_prof_cas(nlev_cas),q2_prof_cas(nlev_cas)
    659687! local:
     
    739767       sens_prof_cas = sens_cas(it_cas2)                                 &
    740768     &          -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
    741        tke_prof_cas = tke_cas(it_cas2)                                   &
    742      &          -frac*(tke_cas(it_cas2)-tke_cas(it_cas1))
     769       tkes_prof_cas = tkes_cas(it_cas2)                                   &
     770     &          -frac*(tkes_cas(it_cas2)-tkes_cas(it_cas1))
    743771       ts_prof_cas = ts_cas(it_cas2)                                     &
    744772     &          -frac*(ts_cas(it_cas2)-ts_cas(it_cas1))
     
    786814        omega_prof_cas(k) = omega_cas(k,it_cas2)                         &
    787815     &          -frac*(omega_cas(k,it_cas2)-omega_cas(k,it_cas1))
     816        tke_prof_cas(k) = tke_cas(k,it_cas2)                         &
     817     &          -frac*(tke_cas(k,it_cas2)-tke_cas(k,it_cas1))
    788818        du_prof_cas(k) = du_cas(k,it_cas2)                               &
    789819     &          -frac*(du_cas(k,it_cas2)-du_cas(k,it_cas1))
     
    833863!**********************************************************************************************
    834864!=====================================================================
    835        SUBROUTINE interp2_case_vertical_std(play,nlev_cas,plev_prof_cas                                &
     865       SUBROUTINE interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                           &
    836866     &         ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas                                       &
    837867     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                              &
    838868     &         ,ug_prof_cas,vg_prof_cas                                                                &
    839869     &         ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                    &
    840      &         ,vitw_prof_cas,omega_prof_cas                                                           &
     870     &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                              &
    841871     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                &
    842872     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
     
    847877     &         ,ug_mod_cas,vg_mod_cas                                                                  &
    848878     &         ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                        &
    849      &         ,w_mod_cas,omega_mod_cas                                                                &
     879     &         ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                    &
    850880     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                      &
    851881     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas        &
     
    870900!       real hq_prof(nlevmax),vq_prof(nlevmax)
    871901 
    872        real play(llm), plev_prof_cas(nlev_cas)
     902       real play(llm), plev(llm+1), plev_prof_cas(nlev_cas)
    873903       real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
    874904       real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
    875905       real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
    876        real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas)
     906       real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas),tke_prof_cas(nlev_cas)
    877907       real temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas)
    878908       real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
     
    887917       real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
    888918       real u_mod_cas(llm),v_mod_cas(llm)
    889        real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm)
     919       real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm),tke_mod_cas(llm+1)
    890920       real temp_nudg_mod_cas(llm),qv_nudg_mod_cas(llm)
    891921       real u_nudg_mod_cas(llm),v_nudg_mod_cas(llm)
     
    899929       real frac,frac1,frac2,fact
    900930 
    901 !       do l = 1, llm
    902 !       print *,'debut interp2, play=',l,play(l)
    903 !       enddo
    904 !      do l = 1, nlev_cas
    905 !      print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l)
    906 !      enddo
     931
     932
     933! for variables defined at the middle of layers
    907934
    908935       do l = 1, llm
     
    932959         endif
    933960
     961
     962
    934963         frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
     964         
    935965         t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
    936966         theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
     
    10291059         ug_mod_cas(l)= ug_prof_cas(nlev_cas)                          !jyg
    10301060         vg_mod_cas(l)= vg_prof_cas(nlev_cas)                          !jyg
    1031          temp_nudg_mod_cas(l)= temp_nudg_prof_cas(nlev_cas)                          !jyg
    1032          qv_nudg_mod_cas(l)= qv_nudg_prof_cas(nlev_cas)                          !jyg
    1033          u_nudg_mod_cas(l)= u_nudg_prof_cas(nlev_cas)                          !jyg
    1034          v_nudg_mod_cas(l)= v_nudg_prof_cas(nlev_cas)                          !jyg
     1061         temp_nudg_mod_cas(l)= temp_nudg_prof_cas(nlev_cas)            !jyg
     1062         qv_nudg_mod_cas(l)= qv_nudg_prof_cas(nlev_cas)                !jyg
     1063         u_nudg_mod_cas(l)= u_nudg_prof_cas(nlev_cas)                  !jyg
     1064         v_nudg_mod_cas(l)= v_nudg_prof_cas(nlev_cas)                  !jyg
    10351065         thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
    10361066         w_mod_cas(l)= 0.0                                             !jyg
     
    10571087       enddo ! l
    10581088
     1089! for variables defined at layer interfaces (EV):
     1090
     1091
     1092       do l = 1, llm+1
     1093
     1094        if (plev(l).ge.plev_prof_cas(nlev_cas)) then
     1095
     1096         mxcalc=l
     1097         k1=0
     1098         k2=0
     1099
     1100         if (plev(l).le.plev_prof_cas(1)) then
     1101
     1102         do k = 1, nlev_cas-1
     1103          if (plev(l).le.plev_prof_cas(k).and. plev(l).gt.plev_prof_cas(k+1)) then
     1104            k1=k
     1105            k2=k+1
     1106          endif
     1107         enddo
     1108
     1109         if (k1.eq.0 .or. k2.eq.0) then
     1110          write(*,*) 'PB! k1, k2 = ',k1,k2
     1111          write(*,*) 'l,plev(l) = ',l,plev(l)/100
     1112         do k = 1, nlev_cas-1
     1113          write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
     1114         enddo
     1115         endif
     1116
     1117         frac = (plev_prof_cas(k2)-plev(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
     1118         tke_mod_cas(l)= tke_prof_cas(k2) - frac*(tke_prof_cas(k2)-tke_prof_cas(k1))
     1119         else !play>plev_prof_cas(1)
     1120         k1=1
     1121         k2=2
     1122         tke_mod_cas(l)= frac1*tke_prof_cas(k1) - frac2*tke_prof_cas(k2)
     1123
     1124         endif ! plev.le.plev_prof_cas(1)
     1125
     1126        else ! above max altitude of forcing file
     1127
     1128         tke_mod_cas(l)=0.0
     1129
     1130        endif ! plev
     1131
     1132       enddo ! l
     1133
     1134
     1135
    10591136          return
    1060           end
     1137          end SUBROUTINE interp2_case_vertical_std
    10611138!*****************************************************************************
    10621139
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_decl_cases.h

    r3605 r3798  
    3737        real th_mod(llm)
    3838
    39         real ts_cur
    40         common /sst_forcing/ts_cur ! also in read_tsurf1d.F
     39        !real ts_cur
     40        !common /sst_forcing/ts_cur ! also in read_tsurf1d.F
    4141!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4242! Declarations specifiques au cas RICO
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_interp_cases.h

    r3605 r3798  
    6262     &             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof    &
    6363     &             ,ht_prof,vt_prof,hq_prof,vq_prof)
    64 
    65         if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
     64! EV: tg instead of ts_cur
     65        if (type_ts_forcing.eq.1) tg = ts_prof !
    6666
    6767! vertical interpolation:
     
    113113!     print *,'llm l omega_profd',llm,l,omega_profd(l)
    114114!     enddo
    115 
    116         if (type_ts_forcing.eq.1) ts_cur = tg_prof ! SST used in read_tsurf1d
     115! EV tg instead of ts_cur
     116        if (type_ts_forcing.eq.1) tg = tg_prof ! SST used
    117117
    118118! vertical interpolation:
     
    206206     &             ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4                            &
    207207     &             ,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
     208!EV tg instead of ts_cur
     209        if (type_ts_forcing.eq.1) tg = tg_prof ! SST used
    210210
    211211! vertical interpolation:
     
    499499     &             ,nlev_sandu                                              &
    500500     &             ,ts_sandu,ts_prof)
    501 
    502         if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
     501! EV tg instead of ts_cur
     502        if (type_ts_forcing.eq.1) tg = ts_prof ! SST used in read_tsurf1d
    503503
    504504! vertical interpolation:
     
    582582     &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof    &
    583583     &             ,ufa_prof,vfa_prof)
    584 
    585         if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
    586 
     584! EV tg instead of ts_cur
     585        if (type_ts_forcing.eq.1) tg = ts_prof ! SST used
    587586! vertical interpolation:
    588587      CALL interp_astex_vertical(play,nlev_astex,plev_profa                 &
     
    675674     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas               &
    676675     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas)
    677 
    678              ts_cur = ts_prof_cas
     676! EV tg instead of ts_cur
     677
     678             tg = ts_prof_cas
    679679             psurf=plev_prof_cas(1)
    680680
     
    850850     &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
    851851     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
    852 
    853              ts_cur = ts_prof_cas
     852! EV tg instead of ts_cur
     853
     854             tg = ts_prof_cas
    854855!            psurf=plev_prof_cas(1)
    855856             psurf=ps_prof_cas
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_1D_read_forc_cases.h

    r3605 r3798  
    875875
    876876! initial and boundary conditions :
    877 !      tsurf = ts_prof_cas
    878       ts_cur = ts_prof_cas
     877!     tsurf = ts_prof_cas
     878! EV tg instead of ts_cur
     879      tg= ts_prof_cas
    879880      psurf=plev_prof_cas(1)
    880881      write(*,*) 'SST initiale: ',tsurf
     
    965966! initial and boundary conditions :
    966967!      tsurf = ts_prof_cas
    967       ts_cur = ts_prof_cas
     968! EV tg instead of ts_cur
     969      tg = ts_prof_cas
    968970      psurf=plev_prof_cas(1)
    969971      write(*,*) 'SST initiale: ',tsurf
     
    10151017      if (forcing_SCM) then
    10161018
    1017          write(*,*),'avant call read_SCM'
    1018          call read_SCM_cas
     1019         write(*,*),'avant call old_read_SCM_cas'
     1020         call old_read_SCM_cas
    10191021         write(*,*) 'Forcing read'
    10201022
     
    10631065! initial and boundary conditions :
    10641066!      tsurf = ts_prof_cas
    1065       ts_cur = ts_prof_cas
     1067! EV tg instead of ts_cur
     1068
     1069      tg = ts_prof_cas
    10661070      psurf=plev_prof_cas(1)
    10671071      write(*,*) 'SST initiale: ',tsurf
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/old_lmdz1d.F90

    r3605 r3798  
    632632!          (phys_state_var_init is called again in physiq)
    633633      read_climoz = 0
    634 !
     634      nsw=6          ! EV et LF: sinon, falb_dir et falb_dif ne peuvent etre alloues
     635
     636
    635637      call phys_state_var_init(read_climoz)
    636638
     
    728730
    729731!Al1 pour SST forced, appell?? depuis ocean_forced_noice
    730       ts_cur = tsurf ! SST used in read_tsurf1d
     732! EV tg instead of ts_cur
     733
     734      tg = tsurf ! SST used in read_tsurf1d
    731735!=====================================================================
    732736! Initialisation de la physique :
     
    791795
    792796        fder=0.
     797        print *, 'snsrf', snsrf
    793798        snsrf(1,:)=snowmass ! masse de neige des sous surface
    794799        qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
     
    841846     end if
    842847
    843 
    844848        print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
    845849     &        ,pctsrf(1,is_oce),pctsrf(1,is_ter)
     
    848852        zpic = zpicinp
    849853        ftsol=tsurf
    850         nsw=6 ! on met le nb de bandes SW=6, pour initialiser
    851               ! 6 albedo, mais on peut quand meme tourner avec
    852               ! moins. Seules les 2 ou 4 premiers seront lus
    853854        falb_dir=albedo
    854855        falb_dif=albedo
     
    913914        v_ancien(1,:)=v(:)
    914915
    915 u10m=0.
    916 v10m=0.
    917 ale_wake=0.
    918 ale_bl_stat=0.
     916        u10m=0.
     917        v10m=0.
     918        ale_wake=0.
     919        ale_bl_stat=0.
    919920 
    920921!------------------------------------------------------------------------
  • LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/scm.F90

    r3605 r3798  
    7575      real :: zcufi    = 1.
    7676      real :: zcvfi    = 1.
    77 
    78 !-      real :: nat_surf
    79 !-      logical :: ok_flux_surf
    80 !-      real :: fsens
    81 !-      real :: flat
    82 !-      real :: tsurf
    83 !-      real :: rugos
    84 !-      real :: qsol(1:2)
    85 !-      real :: qsurf
    86 !-      real :: psurf
    87 !-      real :: zsurf
    88 !-      real :: albedo
    89 !-
    90 !-      real :: time     = 0.
    91 !-      real :: time_ini
    92 !-      real :: xlat
    93 !-      real :: xlon
    94 !-      real :: wtsurf
    95 !-      real :: wqsurf
    96 !-      real :: restart_runoff
    97 !-      real :: xagesno
    98 !-      real :: qsolinp
    99 !-      real :: zpicinp
    100 !-
    10177      real :: fnday
    10278      real :: day, daytime
     
    141117        logical :: forcing_case2   = .false.
    142118        logical :: forcing_SCM   = .false.
    143         integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
    144 !                                                            (cf read_tsurf1d.F)
    145119
    146120!flag forcings
     
    148122        logical :: nudge_thermo=.false.
    149123        logical :: cptadvw=.true.
     124
     125
    150126!=====================================================================
    151127! DECLARATIONS FOR EACH CASE
     
    190166      real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
    191167      real :: d_u_nudge(llm),d_v_nudge(llm)
    192       real :: du_adv(llm),dv_adv(llm)
    193       real :: du_age(llm),dv_age(llm)
     168!      real :: d_u_adv(llm),d_v_adv(llm)
     169      real :: d_u_age(llm),d_v_age(llm)
    194170      real :: alpha
    195171      real :: ttt
     
    248224!
    249225      integer :: it_end ! iteration number of the last call
    250 !Al1
     226!Al1,plev,play,phi,phis,presnivs,
    251227      integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
    252228      data ecrit_slab_oc/-1/
     
    273249      d_u_nudge(:)=0.
    274250      d_v_nudge(:)=0.
    275       du_adv(:)=0.
    276       dv_adv(:)=0.
    277       du_age(:)=0.
    278       dv_age(:)=0.
     251      d_u_adv(:)=0.
     252      d_v_adv(:)=0.
     253      d_u_age(:)=0.
     254      d_v_age(:)=0.
     255     
    279256     
    280257! Initialization of Common turb_forcing
     
    290267! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
    291268!---------------------------------------------------------------------
    292 !Al1
    293269        call conf_unicol
    294270!Al1 moves this gcssold var from common fcg_gcssold to
     
    296272! --------------------------------------------------------------------
    297273        close(1)
    298 !Al1
    299274        write(*,*) 'lmdz1d.def lu => unicol.def'
    300275
     
    302277       year_ini_cas=1997
    303278       ! It is possible that those parameters are run twice.
    304 
    305279       ! A REVOIR : LIRE PEUT ETRE AN MOIS JOUR DIRECETEMENT
     280
     281
    306282       call getin('anneeref',year_ini_cas)
    307283       call getin('dayref',day_deb)
     
    309285       call getin('time_ini',heure_ini_cas)
    310286
    311         type_ts_forcing = 0
    312         IF (nat_surf==0) type_ts_forcing=1 ! SST forcee sur OCEAN
    313         print*,'NATURE DE LA SURFACE ',nat_surf
     287      print*,'NATURE DE LA SURFACE ',nat_surf
    314288!
    315289! Initialization of the logical switch for nudging
     290
    316291     jcode = iflag_nudge
    317292     do i = 1,nudge_max
     
    319294       jcode = jcode/10
    320295     enddo
    321 !---------------------------------------------------------------------
     296!-----------------------------------------------------------------------
    322297!  Definition of the run
    323 !---------------------------------------------------------------------
     298!-----------------------------------------------------------------------
    324299
    325300      call conf_gcm( 99, .TRUE. )
     
    343318      allocate( phy_flic(year_len)) ! Fraction de glace
    344319      phy_flic(:)=0.0
     320
     321
    345322!-----------------------------------------------------------------------
    346323!   Choix du calendrier
     
    373350!      Le numero du jour est dans "day". L heure est traitee separement.
    374351!      La date complete est dans "daytime" (l'unite est le jour).
     352
     353
    375354      if (nday>0) then
    376355         fnday=nday
     
    409388! Initialization of dimensions, geometry and initial state
    410389!---------------------------------------------------------------------
    411 !      call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
     390!     call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
    412391!     but we still need to initialize dimphy module (klon,klev,etc.)  here.
    413392      call init_dimphy1D(1,llm)
     
    433412!          (phys_state_var_init is called again in physiq)
    434413      read_climoz = 0
    435 !
     414      nsw=6
     415
    436416      call phys_state_var_init(read_climoz)
    437417
     
    446426!!! Feedback forcing values for Gateaux differentiation (al1)
    447427!!!=====================================================================
    448 !!! Surface Planck forcing bracketing call radiation
    449 !!      surf_Planck = 0.
    450 !!      surf_Conv   = 0.
    451 !!      write(*,*) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
    452 !!! a mettre dans le lmdz1d.def ou autre
    453 !!
    454428!!
    455429      qsol = qsolinp
     
    469443      ENDIF
    470444      print*,'Flux sol ',fsens,flat
    471 !!      ok_flux_surf=.false.
    472 !!      fsens=-wtsurf*rcpd*rho(1)
    473 !!      flat=-wqsurf*rlvtt*rho(1)
    474 !!!!
    475445
    476446! Vertical discretization and pressure levels at half and mid levels:
     
    496466      plev =ap+bp*psurf
    497467      play = 0.5*(plev(1:llm)+plev(2:llm+1))
    498       zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles
     468      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles.
    499469
    500470      IF (forcing_type .eq. 59) THEN
     
    527497      print*,'mxcalc=',mxcalc
    528498!     print*,'zlay=',zlay(mxcalc)
    529       print*,'play=',play(mxcalc)
    530 
    531 !Al1 pour SST forced, appell?? depuis ocean_forced_noice
    532       ts_cur = tsurf ! SST used in read_tsurf1d
     499!      print*,'play=',play(mxcalc)
     500
     501!! When surface temperature is forced
     502      tg= tsurf ! surface T used in read_tsurf1d
     503
     504
    533505!=====================================================================
    534506! Initialisation de la physique :
     
    546518! airefi,zcufi,zcvfi initialises au debut de ce programme
    547519! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
     520
     521
    548522      day_step = float(nsplit_phys)*day_step/float(iphysiq)
    549523      write (*,*) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
     
    563537     ! e.g. for cell boundaries, which are meaningless in 1D; so pad these
    564538     ! with '0.' when necessary
     539
    565540      call iniphysiq(iim,jjm,llm, &
    566541           1,comm_lmdz, &
     
    650625        zpic = zpicinp
    651626        ftsol=tsurf
    652         nsw=6 ! on met le nb de bandes SW=6, pour initialiser
    653               ! 6 albedo, mais on peut quand meme tourner avec
    654               ! moins. Seules les 2 ou 4 premiers seront lus
    655627        falb_dir=albedo
    656628        falb_dif=albedo
     
    664636        prw_ancien = 0.
    665637!jyg<
    666 !!        pbl_tke(:,:,:)=1.e-8
    667         pbl_tke(:,:,:)=0.
    668         pbl_tke(:,2,:)=1.e-2
    669         PRINT *, ' pbl_tke dans lmdz1d '
    670         if (prt_level .ge. 5) then
    671          DO nsrf = 1,4
    672            PRINT *,'pbl_tke(1,:,',nsrf,') ',pbl_tke(1,:,nsrf)
    673          ENDDO
    674         end if
    675 
     638! Etienne: comment those lines since now the TKE is inialized in 1D_read_forc_cases
     639!!      pbl_tke(:,:,:)=1.e-8
     640!        pbl_tke(:,:,:)=0.
     641!        pbl_tke(:,2,:)=1.e-2
    676642!>jyg
    677 
    678643        rain_fall=0.
    679644        snow_fall=0.
     
    715680        v_ancien(1,:)=v(:)
    716681 
    717 u10m=0.
    718 v10m=0.
    719 ale_wake=0.
    720 ale_bl_stat=0.
     682        u10m=0.
     683        v10m=0.
     684        ale_wake=0.
     685        ale_bl_stat=0.
    721686
    722687!------------------------------------------------------------------------
     
    738703! to be set at some arbitratry convenient values.
    739704!------------------------------------------------------------------------
    740 !Al1 =============== restart option ==========================
     705!Al1 =============== restart option ======================================
    741706        if (.not.restart) then
    742707          iflag_pbl = 5
     
    803768       print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
    804769       print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
    805        print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
     770       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis(1)
    806771! raz for safety
    807772       do l=1,llm
     
    809774       enddo
    810775      endif
    811 !Al1 ================  end restart =================================
     776!======================  end restart =================================
    812777      IF (ecrit_slab_oc.eq.1) then
    813778         open(97,file='div_slab.dat',STATUS='UNKNOWN')
     
    820785       CALL iophys_ini
    821786#endif
     787
     788!=====================================================================
    822789! START OF THE TEMPORAL LOOP :
    823790!=====================================================================
    824791           
    825792      it_end = nint(fnday*day_step)
    826 !test JLD     it_end = 10
    827793      do while(it.le.it_end)
    828794
     
    832798         print*,'PAS DE TEMPS ',timestep
    833799       endif
    834 !Al1 demande de restartphy.nc
    835800       if (it.eq.it_end) lastcall=.True.
    836801
     
    840805
    841806#include "1D_interp_cases.h"
    842    ! Vertical advection
    843 !  call lstendH(llm,nqtot,omega,d_t_vert_adv,d_q_vert_adv,q,temp,u,v,play)
    844 !  print*,'B d_t_adv ',d_t_adv(1:20)*86400
    845 !  print*,'B d_t_vert_adv ',d_t_vert_adv(1:20)*86400
    846 !  print*,'B dt omega ',omega
    847807
    848808!---------------------------------------------------------------------
    849809!  Geopotential :
    850810!---------------------------------------------------------------------
    851 
     811!        phis(1)=zsurf*RG
     812!        phi(1)=phis(1)+RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
    852813        phi(1)=RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
     814
    853815        do l = 1, llm-1
    854816          phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*                           &
    855817     &    (play(l)-play(l+1))/(play(l)+play(l+1))
    856818        enddo
     819
    857820
    858821!---------------------------------------------------------------------
     
    872835      teta=temp*(pzero/play)**rkappa
    873836      do l=2,llm-1
     837        ! vertical tendencies computed as d X / d t = -W  d X / d z
    874838        d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
    875839        d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
    876         d_t_vert_adv(l)=-(w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)))/(pzero/play(l))**rkappa
     840        ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
     841        d_t_vert_adv(l)=-w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))**rkappa
    877842        d_q_vert_adv(l,1)=-w_adv(l)*(q(l+1,1)-q(l-1,1))/(z_adv(l+1)-z_adv(l-1))
    878843      enddo
     844      d_u_adv(:)=d_u_adv(:)+d_u_vert_adv(:)
     845      d_v_adv(:)=d_v_adv(:)+d_v_vert_adv(:)
    879846      d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:)
    880847      d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1)
     
    938905       fcoriolis=2.*sin(rpi*xlat/180.)*romega
    939906
    940 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    941 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    942 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    943 !       if (forcing_radconv .or. forcing_fire) then
    944 !         fcoriolis=0.0
    945 !         dt_cooling=0.0
    946 !         d_t_adv=0.0
    947 !         d_q_adv=0.0
    948 !       endif
    949 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    950 
    951 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    952 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    953 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    954 !      if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice            &
    955 !    &    .or.forcing_amma .or. forcing_type.eq.101) then
    956 !        fcoriolis=0.0 ; ug=0. ; vg=0.
    957 !      endif
    958 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    959 
    960 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    961 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    962 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    963 !      if(forcing_rico) then
    964 !         dt_cooling=0.
    965 !      endif
    966 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    967 
    968 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    969 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    970 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    971 !CRio:Attention modif sp??cifique cas de Caroline
    972 !     if (forcing_type==-1) then
    973 !        fcoriolis=0.
    974 !       
    975 !on calcule dt_cooling
    976 !       do l=1,llm
    977 !       if (play(l).ge.20000.) then
    978 !           dt_cooling(l)=-1.5/86400.
    979 !       elseif ((play(l).ge.10000.).and.((play(l).lt.20000.))) then
    980 !           dt_cooling(l)=-1.5/86400.*(play(l)-10000.)/(10000.)-1./86400.*(20000.-play(l))/10000.*(temp(l)-200.)
    981 !       else
    982 !           dt_cooling(l)=-1.*(temp(l)-200.)/86400.
    983 !       endif
    984 !       enddo
    985 !
    986 !     endif     
    987 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    988 
    989 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    990 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    991 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    992 !     if (forcing_sandu) then
    993 !        ug(1:llm)=u_mod(1:llm)
    994 !        vg(1:llm)=v_mod(1:llm)
    995 !     endif
    996 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    997 
    998907      IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
    999908                                   fcoriolis, xlat,mxcalc
    1000909
    1001 !       print *,'u-ug=',u-ug
    1002 
    1003 !!!!!!!!!!!!!!!!!!!!!!!!
    1004 ! Geostrophic wind
    1005 ! Le calcul ci dessous est insuffisamment precis
    1006 !      du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
    1007 !      dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
    1008 !!!!!!!!!!!!!!!!!!!!!!!!
     910!---------------------------------------------------------------------
     911! Geostrophic forcing
     912!---------------------------------------------------------------------
     913
     914      IF ( forc_geo == 0 ) THEN
     915              d_u_age(1:mxcalc)=0.
     916              d_v_age(1:mxcalc)=0.
     917      ELSE
    1009918       sfdt = sin(0.5*fcoriolis*timestep)
    1010919       cfdt = cos(0.5*fcoriolis*timestep)
    1011 !       print *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep
    1012 !
    1013         du_age(1:mxcalc)= -2.*sfdt/timestep*                                &
     920
     921        d_u_age(1:mxcalc)= -2.*sfdt/timestep*                                &
    1014922     &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
    1015923     &           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
    1016924!!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
    1017925!
    1018        dv_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
     926       d_v_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
    1019927     &          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +                           &
    1020928     &           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
    1021929!!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
    1022 !
    1023 !!!!!!!!!!!!!!!!!!!!!!!!
     930      ENDIF
     931!
     932!---------------------------------------------------------------------
    1024933!  Nudging
    1025 !!!!!!!!!!!!!!!!!!!!!!!!
     934!---------------------------------------------------------------------
    1026935      d_t_nudge(:) = 0.
    1027936      d_u_nudge(:) = 0.
     
    1039948      ENDDO
    1040949
     950!---------------------------------------------------------------------
     951!  Optional outputs
     952!---------------------------------------------------------------------
    1041953#ifdef OUTPUT_PHYS_SCM
    1042954      CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv)
     
    1056968#endif
    1057969
    1058 !
    1059 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1060 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    1061 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1062 !       if (forcing_fire) THEN
    1063 !               print*,'Enlever cette section rapidement'
    1064 !               stop
    1065 !               
    1066 !
    1067 !!let ww=if ( alt le 1100 ) then alt*-0.00001 else 0
    1068 !!let wt=if ( alt le 1100 ) then min( -3.75e-5 , -7.5e-8*alt)  else 0
    1069 !!let wq=if ( alt le 1100 ) then max( 1.5e-8 , 3e-11*alt)  else 0
    1070 !           d_t_adv=0.
    1071 !           d_q_adv=0.
    1072 !           teta=temp*(pzero/play)**rkappa
    1073 !           d_t_adv=0.
    1074 !           d_q_adv=0.
    1075 !           do l=2,llm-1
    1076 !              if (zlay(l)<=1100) then
    1077 !                  wwww=-0.00001*zlay(l)
    1078 !                  d_t_adv(l)=-wwww*(teta(l)-teta(l+1))/(zlay(l)-zlay(l+1)) /(pzero/play(l))**rkappa
    1079 !                  d_q_adv(l,1:2)=-wwww*(q(l,1:2)-q(l+1,1:2))/(zlay(l)-zlay(l+1))
    1080 !                  d_t_adv(l)=d_t_adv(l)+min(-3.75e-5 , -7.5e-8*zlay(l))
    1081 !                  d_q_adv(l,1)=d_q_adv(l,1)+max( 1.5e-8 , 3e-11*zlay(l))
    1082 !              endif
    1083 !           enddo
    1084 !
    1085 !        endif
    1086 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1087 
    1088 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1089 !         call  writefield_phy('dv_age' ,dv_age,llm)
    1090 !         call  writefield_phy('du_age' ,du_age,llm)
    1091 !         call  writefield_phy('du_phys' ,du_phys,llm)
    1092 !         call  writefield_phy('u_tend' ,u,llm)
    1093 !         call  writefield_phy('u_g' ,ug,llm)
    1094 !
    1095 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1096 !! Increment state variables
    1097 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1098970    IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
    1099971
    1100 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1101 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    1102 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1103 ! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h
    1104 ! au dessus de 700hpa, on relaxe vers les profils initiaux
    1105 !     if (forcing_sandu .OR. forcing_astex) then
    1106 !#include "1D_nudge_sandu_astex.h"
    1107 !      else
    1108 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1109972        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
    1110973     &              du_phys(1:mxcalc)                                       &
    1111      &             +du_age(1:mxcalc)+du_adv(1:mxcalc)                       &
     974     &             +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc)                       &
    1112975     &             +d_u_nudge(1:mxcalc) )           
    1113976        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
    1114977     &              dv_phys(1:mxcalc)                                       &
    1115      &             +dv_age(1:mxcalc)+dv_adv(1:mxcalc)                       &
     978     &             +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc)                       &
    1116979     &             +d_v_nudge(1:mxcalc) )
    1117980        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
     
    1125988     &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
    1126989           print* ,'dv_phys=',dv_phys
    1127            print* ,'dv_age=',dv_age
    1128            print* ,'dv_adv=',dv_adv
     990           print* ,'d_v_age=',d_v_age
     991           print* ,'d_v_adv=',d_v_adv
    1129992           print* ,'d_v_nudge=',d_v_nudge
    1130993           print*, v
     
    1134997        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
    1135998     &              dt_phys(1:mxcalc)                                       &
    1136      &             +d_t_adv(1:mxcalc)                                      &
    1137      &             +d_t_nudge(1:mxcalc)                                      &
     999     &             +d_t_adv(1:mxcalc)                                       &
     1000     &             +d_t_nudge(1:mxcalc)                                     &
    11381001     &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
    11391002
    11401003
    1141 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1004!=======================================================================
    11421005!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    1143 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1144 !     endif  ! forcing_sandu or forcing_astex
    1145 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1006!=======================================================================
    11461007
    11471008        teta=temp*(pzero/play)**rkappa
    1148 !
     1009
    11491010!---------------------------------------------------------------------
    11501011!   Nudge soil temperature if requested
     
    11841045
    11851046!  incremente day time
    1186 !        print*,'daytime bef',daytime,1./day_step
    11871047        daytime = daytime+1./day_step
    1188 !Al1dbg
    11891048        day = int(daytime+0.1/day_step)
    11901049!        time = max(daytime-day,0.0)
     
    11921051!cc        time = real(mod(it,day_step))/day_step
    11931052        time = time_ini/24.+real(mod(it,day_step))/day_step
    1194 !        print*,'daytime nxt time',daytime,time
    11951053        it=it+1
    11961054
    11971055      enddo
    11981056
    1199 !Al1
    12001057      if (ecrit_slab_oc.ne.-1) close(97)
    12011058
    12021059!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
    1203 ! -------------------------------------
     1060! ---------------------------------------------------------------------------
    12041061       call dyn1dredem("restart1dyn.nc",                                    &
    12051062     &              plev,play,phi,phis,presnivs,                            &
Note: See TracChangeset for help on using the changeset viewer.