Ignore:
Timestamp:
Mar 29, 2023, 3:14:27 PM (18 months ago)
Author:
lguez
Message:

Sync latest trunk changes to branch LMDZ_ECRad

Location:
LMDZ6/branches/LMDZ_ECRad
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/LMDZ_ECRad

  • LMDZ6/branches/LMDZ_ECRad/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90

    r4104 r4482  
    55
    66!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    7 !Declarations specifiques au cas standard
    8         character*80 :: fich_cas
    9 ! Discr?tisation
    10         integer nlev_cas, nt_cas
    11 
    12 
    13 !profils environnementaux
    14         real, allocatable::  plev_cas(:,:),plevh_cas(:)
    15         real, allocatable::  ap_cas(:),bp_cas(:)
    16 
    17         real, allocatable::  z_cas(:,:),zh_cas(:)
    18         real, allocatable::  t_cas(:,:),q_cas(:,:),qv_cas(:,:),ql_cas(:,:),qi_cas(:,:),rh_cas(:,:)
    19         real, allocatable::  th_cas(:,:),thv_cas(:,:),thl_cas(:,:),rv_cas(:,:)
    20         real, allocatable::  u_cas(:,:),v_cas(:,:),vitw_cas(:,:),omega_cas(:,:),tke_cas(:,:)
    21 
    22 !forcing
    23         real, allocatable::  ht_cas(:,:),vt_cas(:,:),dt_cas(:,:),dtrad_cas(:,:)
    24         real, allocatable::  hth_cas(:,:),vth_cas(:,:),dth_cas(:,:)
    25         real, allocatable::  hq_cas(:,:),vq_cas(:,:),dq_cas(:,:)
    26         real, allocatable::  hr_cas(:,:),vr_cas(:,:),dr_cas(:,:)
    27         real, allocatable::  hu_cas(:,:),vu_cas(:,:),du_cas(:,:)
    28         real, allocatable::  hv_cas(:,:),vv_cas(:,:),dv_cas(:,:)
    29         real, allocatable::  ug_cas(:,:),vg_cas(:,:)
    30         real, allocatable::  temp_nudg_cas(:,:),qv_nudg_cas(:,:),u_nudg_cas(:,:),v_nudg_cas(:,:)
    31         real, allocatable::  invtau_temp_nudg_cas(:,:),invtau_qv_nudg_cas(:,:),invtau_u_nudg_cas(:,:),invtau_v_nudg_cas(:,:)
    32         real, allocatable::  lat_cas(:),sens_cas(:),ts_cas(:),ps_cas(:),ustar_cas(:)
    33         real, allocatable::  uw_cas(:,:),vw_cas(:,:),q1_cas(:,:),q2_cas(:,:),tkes_cas(:)
    34 
    35 !champs interpoles
    36         real, allocatable::  plev_prof_cas(:)
    37         real, allocatable::  t_prof_cas(:)
    38         real, allocatable::  theta_prof_cas(:)
    39         real, allocatable::  thl_prof_cas(:)
    40         real, allocatable::  thv_prof_cas(:)
    41         real, allocatable::  q_prof_cas(:)
    42         real, allocatable::  qv_prof_cas(:)
    43         real, allocatable::  ql_prof_cas(:)
    44         real, allocatable::  qi_prof_cas(:)
    45         real, allocatable::  rh_prof_cas(:)
    46         real, allocatable::  rv_prof_cas(:)
    47         real, allocatable::  u_prof_cas(:)
    48         real, allocatable::  v_prof_cas(:)       
    49         real, allocatable::  vitw_prof_cas(:)
    50         real, allocatable::  omega_prof_cas(:)
    51         real, allocatable::  tke_prof_cas(:)
    52         real, allocatable::  ug_prof_cas(:)
    53         real, allocatable::  vg_prof_cas(:)
    54         real, allocatable::  temp_nudg_prof_cas(:),qv_nudg_prof_cas(:),u_nudg_prof_cas(:),v_nudg_prof_cas(:)
    55         real, allocatable::  invtau_temp_nudg_prof_cas(:),invtau_qv_nudg_prof_cas(:),invtau_u_nudg_prof_cas(:),invtau_v_nudg_prof_cas(:)
    56 
    57         real, allocatable::  ht_prof_cas(:)
    58         real, allocatable::  hth_prof_cas(:)
    59         real, allocatable::  hq_prof_cas(:)
    60         real, allocatable::  vt_prof_cas(:)
    61         real, allocatable::  vth_prof_cas(:)
    62         real, allocatable::  vq_prof_cas(:)
    63         real, allocatable::  dt_prof_cas(:)
    64         real, allocatable::  dth_prof_cas(:)
    65         real, allocatable::  dtrad_prof_cas(:)
    66         real, allocatable::  dq_prof_cas(:)
    67         real, allocatable::  hu_prof_cas(:)
    68         real, allocatable::  hv_prof_cas(:)
    69         real, allocatable::  vu_prof_cas(:)
    70         real, allocatable::  vv_prof_cas(:)
    71         real, allocatable::  du_prof_cas(:)
    72         real, allocatable::  dv_prof_cas(:)
    73         real, allocatable::  uw_prof_cas(:)
    74         real, allocatable::  vw_prof_cas(:)
    75         real, allocatable::  q1_prof_cas(:)
    76         real, allocatable::  q2_prof_cas(:)
    77 
    78 
    79         real o3_cas,lat_prof_cas,sens_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas,tkes_prof_cas
    80         real orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough,rugos_cas,sand_cas,clay_cas
    81      
     7  !Declarations specifiques au cas standard
     8  character*80 :: fich_cas
     9  ! Discr?tisation
     10  integer nlev_cas, nt_cas
     11
     12
     13  !profils environnementaux
     14  real, allocatable::  plev_cas(:,:),plevh_cas(:)
     15  real, allocatable::  ap_cas(:),bp_cas(:)
     16
     17  real, allocatable::  z_cas(:,:),zh_cas(:)
     18  real, allocatable::  t_cas(:,:),q_cas(:,:),qv_cas(:,:),ql_cas(:,:),qi_cas(:,:),rh_cas(:,:)
     19  real, allocatable::  th_cas(:,:),thv_cas(:,:),thl_cas(:,:),rv_cas(:,:)
     20  real, allocatable::  u_cas(:,:),v_cas(:,:),vitw_cas(:,:),omega_cas(:,:),tke_cas(:,:)
     21
     22  !forcing
     23  real, allocatable::  ht_cas(:,:),vt_cas(:,:),dt_cas(:,:),dtrad_cas(:,:)
     24  real, allocatable::  hth_cas(:,:),vth_cas(:,:),dth_cas(:,:)
     25  real, allocatable::  hq_cas(:,:),vq_cas(:,:),dq_cas(:,:)
     26  real, allocatable::  hr_cas(:,:),vr_cas(:,:),dr_cas(:,:)
     27  real, allocatable::  hu_cas(:,:),vu_cas(:,:),du_cas(:,:)
     28  real, allocatable::  hv_cas(:,:),vv_cas(:,:),dv_cas(:,:)
     29  real, allocatable::  ug_cas(:,:),vg_cas(:,:)
     30  real, allocatable::  temp_nudg_cas(:,:),qv_nudg_cas(:,:),u_nudg_cas(:,:),v_nudg_cas(:,:)
     31  real, allocatable::  invtau_temp_nudg_cas(:,:),invtau_qv_nudg_cas(:,:),invtau_u_nudg_cas(:,:),invtau_v_nudg_cas(:,:)
     32  real, allocatable::  lat_cas(:),sens_cas(:),tskin_cas(:),ts_cas(:),ps_cas(:),ustar_cas(:)
     33  real, allocatable::  uw_cas(:,:),vw_cas(:,:),q1_cas(:,:),q2_cas(:,:),tkes_cas(:)
     34
     35  !champs interpoles
     36  real, allocatable::  plev_prof_cas(:)
     37  real, allocatable::  t_prof_cas(:)
     38  real, allocatable::  theta_prof_cas(:)
     39  real, allocatable::  thl_prof_cas(:)
     40  real, allocatable::  thv_prof_cas(:)
     41  real, allocatable::  q_prof_cas(:)
     42  real, allocatable::  qv_prof_cas(:)
     43  real, allocatable::  ql_prof_cas(:)
     44  real, allocatable::  qi_prof_cas(:)
     45  real, allocatable::  rh_prof_cas(:)
     46  real, allocatable::  rv_prof_cas(:)
     47  real, allocatable::  u_prof_cas(:)
     48  real, allocatable::  v_prof_cas(:)       
     49  real, allocatable::  vitw_prof_cas(:)
     50  real, allocatable::  omega_prof_cas(:)
     51  real, allocatable::  tke_prof_cas(:)
     52  real, allocatable::  ug_prof_cas(:)
     53  real, allocatable::  vg_prof_cas(:)
     54  real, allocatable::  temp_nudg_prof_cas(:),qv_nudg_prof_cas(:),u_nudg_prof_cas(:),v_nudg_prof_cas(:)
     55  real, allocatable::  invtau_temp_nudg_prof_cas(:),invtau_qv_nudg_prof_cas(:),invtau_u_nudg_prof_cas(:),invtau_v_nudg_prof_cas(:)
     56
     57  real, allocatable::  ht_prof_cas(:)
     58  real, allocatable::  hth_prof_cas(:)
     59  real, allocatable::  hq_prof_cas(:)
     60  real, allocatable::  vt_prof_cas(:)
     61  real, allocatable::  vth_prof_cas(:)
     62  real, allocatable::  vq_prof_cas(:)
     63  real, allocatable::  dt_prof_cas(:)
     64  real, allocatable::  dth_prof_cas(:)
     65  real, allocatable::  dtrad_prof_cas(:)
     66  real, allocatable::  dq_prof_cas(:)
     67  real, allocatable::  hu_prof_cas(:)
     68  real, allocatable::  hv_prof_cas(:)
     69  real, allocatable::  vu_prof_cas(:)
     70  real, allocatable::  vv_prof_cas(:)
     71  real, allocatable::  du_prof_cas(:)
     72  real, allocatable::  dv_prof_cas(:)
     73  real, allocatable::  uw_prof_cas(:)
     74  real, allocatable::  vw_prof_cas(:)
     75  real, allocatable::  q1_prof_cas(:)
     76  real, allocatable::  q2_prof_cas(:)
     77
     78
     79  real o3_cas,lat_prof_cas,sens_prof_cas,ts_prof_cas,tskin_prof_cas,ps_prof_cas,ustar_prof_cas,tkes_prof_cas
     80  real orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,heat_rough,rugos_cas,sand_cas,clay_cas
     81
    8282
    8383
    … …  
    8585
    8686
    87 !**********************************************************************************************
    88 SUBROUTINE read_SCM_cas
    89       implicit none
     87  !**********************************************************************************************
     88  SUBROUTINE read_SCM_cas
     89    use netcdf, only: nf90_get_var
     90    implicit none
    9091
    9192#include "netcdf.inc"
    9293#include "date_cas.h"
    9394
    94       INTEGER nid,rid,ierr
    95       INTEGER ii,jj,timeid
    96       REAL, ALLOCATABLE :: time_val(:)
    97 
    98       fich_cas='cas.nc'
    99       print*,'fich_cas ',fich_cas
    100       ierr = NF_OPEN(fich_cas,NF_NOWRITE,nid)
    101       print*,'fich_cas,NF_NOWRITE,nid ',fich_cas,NF_NOWRITE,nid
    102       if (ierr.NE.NF_NOERR) then
    103          write(*,*) 'ERROR: GROS Pb opening forcings nc file '
    104          write(*,*) NF_STRERROR(ierr)
    105          stop ""
    106       endif
    107 !.......................................................................
    108       ierr=NF_INQ_DIMID(nid,'lat',rid)
    109       IF (ierr.NE.NF_NOERR) THEN
    110          print*, 'Oh probleme lecture dimension lat'
    111       ENDIF
    112       ierr=NF_INQ_DIMLEN(nid,rid,ii)
    113       print*,'OK1 read2: nid,rid,lat',nid,rid,ii
    114 !.......................................................................
    115       ierr=NF_INQ_DIMID(nid,'lon',rid)
    116       IF (ierr.NE.NF_NOERR) THEN
    117          print*, 'Oh probleme lecture dimension lon'
    118       ENDIF
    119       ierr=NF_INQ_DIMLEN(nid,rid,jj)
    120       print*,'OK2 read2: nid,rid,lat',nid,rid,jj
    121 !.......................................................................
    122       ierr=NF_INQ_DIMID(nid,'lev',rid)
    123       IF (ierr.NE.NF_NOERR) THEN
    124          print*, 'Oh probleme lecture dimension nlev'
    125       ENDIF
    126       ierr=NF_INQ_DIMLEN(nid,rid,nlev_cas)
    127       print*,'OK3 read2: nid,rid,nlev_cas',nid,rid,nlev_cas
    128       IF ( .NOT. ( nlev_cas > 10 .AND. nlev_cas < 200000 )) THEN
    129               print*,'Valeur de nlev_cas peu probable'
    130               STOP
    131       ENDIF
    132 !.......................................................................
    133       ierr=NF_INQ_DIMID(nid,'time',rid)
    134       nt_cas=0
    135       IF (ierr.NE.NF_NOERR) THEN
    136         stop 'Oh probleme lecture dimension time'
    137       ENDIF
    138       ierr=NF_INQ_DIMLEN(nid,rid,nt_cas)
    139       print*,'OK4 read2: nid,rid,nt_cas',nid,rid,nt_cas
    140 ! Lecture de l'axe des temps
    141       print*,'LECTURE DU TEMPS'
    142       ierr=NF_INQ_VARID(nid,'time',timeid)
    143          if(ierr/=NF_NOERR) then
    144            print *,'Variable time manquante dans cas.nc:'
    145            ierr=NF_NOERR
    146          else
    147                  allocate(time_val(nt_cas))
    148 #ifdef NC_DOUBLE
    149          ierr = NF_GET_VAR_DOUBLE(nid,timeid,time_val)
    150 #else
    151            ierr = NF_GET_VAR_REAL(nid,timeid,time_val)
    152 #endif
    153            if(ierr/=NF_NOERR) then
    154               print *,'Pb a la lecture de time cas.nc: '
    155            endif
    156    endif
    157    IF (nt_cas>1) THEN
    158            pdt_cas=time_val(2)-time_val(1)
    159    ELSE
    160            pdt_cas=0.
    161    ENDIF
     95    INTEGER nid,rid,ierr
     96    INTEGER ii,jj,timeid
     97    REAL, ALLOCATABLE :: time_val(:)
     98
     99    fich_cas='cas.nc'
     100    print*,'fich_cas ',fich_cas
     101    ierr = NF_OPEN(fich_cas,NF_NOWRITE,nid)
     102    print*,'fich_cas,NF_NOWRITE,nid ',fich_cas,NF_NOWRITE,nid
     103    if (ierr.NE.NF_NOERR) then
     104       write(*,*) 'ERROR: GROS Pb opening forcings nc file '
     105       write(*,*) NF_STRERROR(ierr)
     106       stop ""
     107    endif
     108    !.......................................................................
     109    ierr=NF_INQ_DIMID(nid,'lat',rid)
     110    IF (ierr.NE.NF_NOERR) THEN
     111       print*, 'Oh probleme lecture dimension lat'
     112    ENDIF
     113    ierr=NF_INQ_DIMLEN(nid,rid,ii)
     114    print*,'OK1 read_SCM_cas: nid,rid,lat',nid,rid,ii
     115    !.......................................................................
     116    ierr=NF_INQ_DIMID(nid,'lon',rid)
     117    IF (ierr.NE.NF_NOERR) THEN
     118       print*, 'Oh probleme lecture dimension lon'
     119    ENDIF
     120    ierr=NF_INQ_DIMLEN(nid,rid,jj)
     121    print*,'OK2 read_SCM_cas: nid,rid,lat',nid,rid,jj
     122    !.......................................................................
     123    ierr=NF_INQ_DIMID(nid,'lev',rid)
     124    IF (ierr.NE.NF_NOERR) THEN
     125       print*, 'Oh probleme lecture dimension nlev'
     126    ENDIF
     127    ierr=NF_INQ_DIMLEN(nid,rid,nlev_cas)
     128    print*,'OK3 read_SCM_cas: nid,rid,nlev_cas',nid,rid,nlev_cas
     129    IF ( .NOT. ( nlev_cas > 10 .AND. nlev_cas < 200000 )) THEN
     130       print*,'Valeur de nlev_cas peu probable'
     131       STOP
     132    ENDIF
     133    !.......................................................................
     134    ierr=NF_INQ_DIMID(nid,'time',rid)
     135    nt_cas=0
     136    IF (ierr.NE.NF_NOERR) THEN
     137       stop 'Oh probleme lecture dimension time'
     138    ENDIF
     139    ierr=NF_INQ_DIMLEN(nid,rid,nt_cas)
     140    print*,'OK4 read_SCM_cas: nid,rid,nt_cas',nid,rid,nt_cas
     141    ! Lecture de l'axe des temps
     142    print*,'LECTURE DU TEMPS'
     143    ierr=NF_INQ_VARID(nid,'time',timeid)
     144    if(ierr/=NF_NOERR) then
     145       print *,'Variable time manquante dans cas.nc:'
     146       ierr=NF_NOERR
     147    else
     148       allocate(time_val(nt_cas))
     149       ierr = NF90_GET_VAR(nid,timeid,time_val)
     150       if(ierr/=NF_NOERR) then
     151          print *,'A Pb a la lecture de time cas.nc: '
     152       endif
     153    endif
     154    IF (nt_cas>1) THEN
     155       pdt_cas=time_val(2)-time_val(1)
     156    ELSE
     157       pdt_cas=0.
     158    ENDIF
    162159
    163160
    164161!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    165 !profils moyens:
    166         allocate(plev_cas(nlev_cas,nt_cas),plevh_cas(nlev_cas+1))       
    167         allocate(z_cas(nlev_cas,nt_cas),zh_cas(nlev_cas+1))
    168         allocate(ap_cas(nlev_cas+1),bp_cas(nt_cas+1))
    169         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), &
    170              qi_cas(nlev_cas,nt_cas),rh_cas(nlev_cas,nt_cas))
    171         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))
    172         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))
    173         allocate(tke_cas(nlev_cas,nt_cas))
    174 !forcing
    175         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))
    176         allocate(hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas),dq_cas(nlev_cas,nt_cas))
    177         allocate(hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas),dth_cas(nlev_cas,nt_cas))
    178         allocate(hr_cas(nlev_cas,nt_cas),vr_cas(nlev_cas,nt_cas),dr_cas(nlev_cas,nt_cas))
    179         allocate(hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas),du_cas(nlev_cas,nt_cas))
    180         allocate(hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas),dv_cas(nlev_cas,nt_cas))
    181         allocate(ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas))
    182         allocate(temp_nudg_cas(nlev_cas,nt_cas),qv_nudg_cas(nlev_cas,nt_cas))
    183         allocate(u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas))
    184         allocate(invtau_temp_nudg_cas(nlev_cas,nt_cas),invtau_qv_nudg_cas(nlev_cas,nt_cas))
    185         allocate(invtau_u_nudg_cas(nlev_cas,nt_cas),invtau_v_nudg_cas(nlev_cas,nt_cas))
    186         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))
    187         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))
    188 
    189 
    190 
    191 !champs interpoles
    192         allocate(plev_prof_cas(nlev_cas))
    193         allocate(t_prof_cas(nlev_cas))
    194         allocate(theta_prof_cas(nlev_cas))
    195         allocate(thl_prof_cas(nlev_cas))
    196         allocate(thv_prof_cas(nlev_cas))
    197         allocate(q_prof_cas(nlev_cas))
    198         allocate(qv_prof_cas(nlev_cas))
    199         allocate(ql_prof_cas(nlev_cas))
    200         allocate(qi_prof_cas(nlev_cas))
    201         allocate(rh_prof_cas(nlev_cas))
    202         allocate(rv_prof_cas(nlev_cas))
    203         allocate(u_prof_cas(nlev_cas))
    204         allocate(v_prof_cas(nlev_cas))
    205         allocate(vitw_prof_cas(nlev_cas))
    206         allocate(omega_prof_cas(nlev_cas))
    207         allocate(tke_prof_cas(nlev_cas))
    208         allocate(ug_prof_cas(nlev_cas))
    209         allocate(vg_prof_cas(nlev_cas))
    210         allocate(temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas))
    211         allocate(u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas))
    212         allocate(invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas))
    213         allocate(invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas))
    214         allocate(ht_prof_cas(nlev_cas))
    215         allocate(hth_prof_cas(nlev_cas))
    216         allocate(hq_prof_cas(nlev_cas))
    217         allocate(hu_prof_cas(nlev_cas))
    218         allocate(hv_prof_cas(nlev_cas))
    219         allocate(vt_prof_cas(nlev_cas))
    220         allocate(vth_prof_cas(nlev_cas))
    221         allocate(vq_prof_cas(nlev_cas))
    222         allocate(vu_prof_cas(nlev_cas))
    223         allocate(vv_prof_cas(nlev_cas))
    224         allocate(dt_prof_cas(nlev_cas))
    225         allocate(dth_prof_cas(nlev_cas))
    226         allocate(dtrad_prof_cas(nlev_cas))
    227         allocate(dq_prof_cas(nlev_cas))
    228         allocate(du_prof_cas(nlev_cas))
    229         allocate(dv_prof_cas(nlev_cas))
    230         allocate(uw_prof_cas(nlev_cas))
    231         allocate(vw_prof_cas(nlev_cas))
    232         allocate(q1_prof_cas(nlev_cas))
    233         allocate(q2_prof_cas(nlev_cas))
    234 
    235         print*,'Allocations OK'
    236         CALL read_SCM (nid,nlev_cas,nt_cas,                                                                     &
    237      &     ap_cas,bp_cas,z_cas,plev_cas,zh_cas,plevh_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,                   &
    238      &     ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,tke_cas,ug_cas,vg_cas,                            &
    239      &     temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas,                                                     &
    240      &     invtau_temp_nudg_cas,invtau_qv_nudg_cas,invtau_u_nudg_cas,invtau_v_nudg_cas,                         &
    241      &     du_cas,hu_cas,vu_cas,                                                                                &
    242      &     dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,               &
    243      &     dr_cas,hr_cas,vr_cas,dtrad_cas,sens_cas,lat_cas,ts_cas,ps_cas,ustar_cas,tkes_cas,                      &
    244      &     uw_cas,vw_cas,q1_cas,q2_cas,orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,heat_rough, &
    245      &     o3_cas,rugos_cas,clay_cas,sand_cas)
    246         print*,'read_SCM cas OK'
    247         do ii=1,nlev_cas
    248         print*,'apres read2_SCM, plev_cas=',ii,plev_cas(ii,1)
    249         !print*,'apres read_SCM, plev_cas=',ii,omega_cas(ii,nt_cas/2+1)
    250         enddo
    251 
    252 
    253 END SUBROUTINE read_SCM_cas
     162    !profils moyens:
     163    allocate(plev_cas(nlev_cas,nt_cas),plevh_cas(nlev_cas+1))       
     164    allocate(z_cas(nlev_cas,nt_cas),zh_cas(nlev_cas+1))
     165    allocate(ap_cas(nlev_cas+1),bp_cas(nt_cas+1))
     166    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), &
     167         qi_cas(nlev_cas,nt_cas),rh_cas(nlev_cas,nt_cas))
     168    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))
     169    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    allocate(tke_cas(nlev_cas,nt_cas))
     171    !forcing
     172    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))
     173    allocate(hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas),dq_cas(nlev_cas,nt_cas))
     174    allocate(hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas),dth_cas(nlev_cas,nt_cas))
     175    allocate(hr_cas(nlev_cas,nt_cas),vr_cas(nlev_cas,nt_cas),dr_cas(nlev_cas,nt_cas))
     176    allocate(hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas),du_cas(nlev_cas,nt_cas))
     177    allocate(hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas),dv_cas(nlev_cas,nt_cas))
     178    allocate(ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas))
     179    allocate(temp_nudg_cas(nlev_cas,nt_cas),qv_nudg_cas(nlev_cas,nt_cas))
     180    allocate(u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas))
     181    allocate(invtau_temp_nudg_cas(nlev_cas,nt_cas),invtau_qv_nudg_cas(nlev_cas,nt_cas))
     182    allocate(invtau_u_nudg_cas(nlev_cas,nt_cas),invtau_v_nudg_cas(nlev_cas,nt_cas))
     183    allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas),tskin_cas(nt_cas),ps_cas(nt_cas),ustar_cas(nt_cas),tkes_cas(nt_cas))
     184    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))
     185
     186
     187
     188    !champs interpoles
     189    allocate(plev_prof_cas(nlev_cas))
     190    allocate(t_prof_cas(nlev_cas))
     191    allocate(theta_prof_cas(nlev_cas))
     192    allocate(thl_prof_cas(nlev_cas))
     193    allocate(thv_prof_cas(nlev_cas))
     194    allocate(q_prof_cas(nlev_cas))
     195    allocate(qv_prof_cas(nlev_cas))
     196    allocate(ql_prof_cas(nlev_cas))
     197    allocate(qi_prof_cas(nlev_cas))
     198    allocate(rh_prof_cas(nlev_cas))
     199    allocate(rv_prof_cas(nlev_cas))
     200    allocate(u_prof_cas(nlev_cas))
     201    allocate(v_prof_cas(nlev_cas))
     202    allocate(vitw_prof_cas(nlev_cas))
     203    allocate(omega_prof_cas(nlev_cas))
     204    allocate(tke_prof_cas(nlev_cas))
     205    allocate(ug_prof_cas(nlev_cas))
     206    allocate(vg_prof_cas(nlev_cas))
     207    allocate(temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas))
     208    allocate(u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas))
     209    allocate(invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas))
     210    allocate(invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas))
     211    allocate(ht_prof_cas(nlev_cas))
     212    allocate(hth_prof_cas(nlev_cas))
     213    allocate(hq_prof_cas(nlev_cas))
     214    allocate(hu_prof_cas(nlev_cas))
     215    allocate(hv_prof_cas(nlev_cas))
     216    allocate(vt_prof_cas(nlev_cas))
     217    allocate(vth_prof_cas(nlev_cas))
     218    allocate(vq_prof_cas(nlev_cas))
     219    allocate(vu_prof_cas(nlev_cas))
     220    allocate(vv_prof_cas(nlev_cas))
     221    allocate(dt_prof_cas(nlev_cas))
     222    allocate(dth_prof_cas(nlev_cas))
     223    allocate(dtrad_prof_cas(nlev_cas))
     224    allocate(dq_prof_cas(nlev_cas))
     225    allocate(du_prof_cas(nlev_cas))
     226    allocate(dv_prof_cas(nlev_cas))
     227    allocate(uw_prof_cas(nlev_cas))
     228    allocate(vw_prof_cas(nlev_cas))
     229    allocate(q1_prof_cas(nlev_cas))
     230    allocate(q2_prof_cas(nlev_cas))
     231
     232    print*,'Allocations OK'
     233    CALL read_SCM (nid,nlev_cas,nt_cas,                                                                     &
     234         ap_cas,bp_cas,z_cas,plev_cas,zh_cas,plevh_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,                   &
     235         ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,tke_cas,ug_cas,vg_cas,                            &
     236         temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas,                                                     &
     237         invtau_temp_nudg_cas,invtau_qv_nudg_cas,invtau_u_nudg_cas,invtau_v_nudg_cas,                         &
     238         du_cas,hu_cas,vu_cas,                                                                                &
     239         dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,               &
     240         dr_cas,hr_cas,vr_cas,dtrad_cas,sens_cas,lat_cas,ts_cas,tskin_cas,ps_cas,ustar_cas,tkes_cas,            &
     241         uw_cas,vw_cas,q1_cas,q2_cas,orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,heat_rough, &
     242         o3_cas,rugos_cas,clay_cas,sand_cas)
     243    print*,'read_SCM cas OK'
     244    do ii=1,nlev_cas
     245       print*,'apres read_SCM_cas, plev_cas=',ii,plev_cas(ii,1)
     246       !print*,'apres read_SCM, plev_cas=',ii,omega_cas(ii,nt_cas/2+1)
     247    enddo
     248
     249
     250  END SUBROUTINE read_SCM_cas
    254251
    255252
    256253!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    257 SUBROUTINE deallocate2_1D_cases
    258 !profils environnementaux:
    259         deallocate(plev_cas,plevh_cas)
    260        
    261         deallocate(z_cas,zh_cas)
    262         deallocate(ap_cas,bp_cas)
    263         deallocate(t_cas,q_cas,qv_cas,ql_cas,qi_cas,rh_cas)
    264         deallocate(th_cas,thl_cas,thv_cas,rv_cas)
    265         deallocate(u_cas,v_cas,vitw_cas,omega_cas,tke_cas)
    266        
    267 !forcing
    268         deallocate(ht_cas,vt_cas,dt_cas,dtrad_cas)
    269         deallocate(hq_cas,vq_cas,dq_cas)
    270         deallocate(hth_cas,vth_cas,dth_cas)
    271         deallocate(hr_cas,vr_cas,dr_cas)
    272         deallocate(hu_cas,vu_cas,du_cas)
    273         deallocate(hv_cas,vv_cas,dv_cas)
    274         deallocate(ug_cas)
    275         deallocate(vg_cas)
    276         deallocate(lat_cas,sens_cas,ts_cas,ps_cas,ustar_cas,tkes_cas,uw_cas,vw_cas,q1_cas,q2_cas)
    277 
    278 !champs interpoles
    279         deallocate(plev_prof_cas)
    280         deallocate(t_prof_cas)
    281         deallocate(theta_prof_cas)
    282         deallocate(thl_prof_cas)
    283         deallocate(thv_prof_cas)
    284         deallocate(q_prof_cas)
    285         deallocate(qv_prof_cas)
    286         deallocate(ql_prof_cas)
    287         deallocate(qi_prof_cas)
    288         deallocate(rh_prof_cas)
    289         deallocate(rv_prof_cas)
    290         deallocate(u_prof_cas)
    291         deallocate(v_prof_cas)
    292         deallocate(vitw_prof_cas)
    293         deallocate(omega_prof_cas)
    294         deallocate(tke_prof_cas)
    295         deallocate(ug_prof_cas)
    296         deallocate(vg_prof_cas)
    297         deallocate(temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas)
    298         deallocate(invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas)
    299         deallocate(ht_prof_cas)
    300         deallocate(hq_prof_cas)
    301         deallocate(hu_prof_cas)
    302         deallocate(hv_prof_cas)
    303         deallocate(vt_prof_cas)
    304         deallocate(vq_prof_cas)
    305         deallocate(vu_prof_cas)
    306         deallocate(vv_prof_cas)
    307         deallocate(dt_prof_cas)
    308         deallocate(dtrad_prof_cas)
    309         deallocate(dq_prof_cas)
    310         deallocate(du_prof_cas)
    311         deallocate(dv_prof_cas)
    312         deallocate(t_prof_cas)
    313         deallocate(u_prof_cas)
    314         deallocate(v_prof_cas)
    315         deallocate(uw_prof_cas)
    316         deallocate(vw_prof_cas)
    317         deallocate(q1_prof_cas)
    318         deallocate(q2_prof_cas)
    319 
    320 END SUBROUTINE deallocate2_1D_cases
    321 
    322 
    323 !=====================================================================
    324       SUBROUTINE read_SCM(nid,nlevel,ntime,                                       &
    325      &     ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,tke,ug,vg,&
    326      &     temp_nudg,qv_nudg,u_nudg,v_nudg,                                        &
    327      &     invtau_temp_nudg,invtau_qv_nudg,invtau_u_nudg,invtau_v_nudg,             &
    328      &     du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq,                                    &
    329      &     dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,ps,ustar,tkes,uw,vw,q1,q2,       &
    330      &     orog_cas,albedo_cas,emiss_cas,t_skin_cas,q_skin_cas,mom_rough,          &
    331      &     heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas)
    332 
    333 !program reading forcing of the case study
    334       implicit none
     254  SUBROUTINE deallocate2_1D_cases
     255    !profils environnementaux:
     256    deallocate(plev_cas,plevh_cas)
     257
     258    deallocate(z_cas,zh_cas)
     259    deallocate(ap_cas,bp_cas)
     260    deallocate(t_cas,q_cas,qv_cas,ql_cas,qi_cas,rh_cas)
     261    deallocate(th_cas,thl_cas,thv_cas,rv_cas)
     262    deallocate(u_cas,v_cas,vitw_cas,omega_cas,tke_cas)
     263
     264    !forcing
     265    deallocate(ht_cas,vt_cas,dt_cas,dtrad_cas)
     266    deallocate(hq_cas,vq_cas,dq_cas)
     267    deallocate(hth_cas,vth_cas,dth_cas)
     268    deallocate(hr_cas,vr_cas,dr_cas)
     269    deallocate(hu_cas,vu_cas,du_cas)
     270    deallocate(hv_cas,vv_cas,dv_cas)
     271    deallocate(ug_cas)
     272    deallocate(vg_cas)
     273    deallocate(lat_cas,sens_cas,tskin_cas,ts_cas,ps_cas,ustar_cas,tkes_cas,uw_cas,vw_cas,q1_cas,q2_cas)
     274
     275    !champs interpoles
     276    deallocate(plev_prof_cas)
     277    deallocate(t_prof_cas)
     278    deallocate(theta_prof_cas)
     279    deallocate(thl_prof_cas)
     280    deallocate(thv_prof_cas)
     281    deallocate(q_prof_cas)
     282    deallocate(qv_prof_cas)
     283    deallocate(ql_prof_cas)
     284    deallocate(qi_prof_cas)
     285    deallocate(rh_prof_cas)
     286    deallocate(rv_prof_cas)
     287    deallocate(u_prof_cas)
     288    deallocate(v_prof_cas)
     289    deallocate(vitw_prof_cas)
     290    deallocate(omega_prof_cas)
     291    deallocate(tke_prof_cas)
     292    deallocate(ug_prof_cas)
     293    deallocate(vg_prof_cas)
     294    deallocate(temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas)
     295    deallocate(invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas)
     296    deallocate(ht_prof_cas)
     297    deallocate(hq_prof_cas)
     298    deallocate(hu_prof_cas)
     299    deallocate(hv_prof_cas)
     300    deallocate(vt_prof_cas)
     301    deallocate(vq_prof_cas)
     302    deallocate(vu_prof_cas)
     303    deallocate(vv_prof_cas)
     304    deallocate(dt_prof_cas)
     305    deallocate(dtrad_prof_cas)
     306    deallocate(dq_prof_cas)
     307    deallocate(du_prof_cas)
     308    deallocate(dv_prof_cas)
     309    deallocate(t_prof_cas)
     310    deallocate(u_prof_cas)
     311    deallocate(v_prof_cas)
     312    deallocate(uw_prof_cas)
     313    deallocate(vw_prof_cas)
     314    deallocate(q1_prof_cas)
     315    deallocate(q2_prof_cas)
     316
     317  END SUBROUTINE deallocate2_1D_cases
     318
     319
     320  !=====================================================================
     321  SUBROUTINE read_SCM(nid,nlevel,ntime,                                       &
     322       ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,tke,ug,vg,&
     323       temp_nudg,qv_nudg,u_nudg,v_nudg,                                        &
     324       invtau_temp_nudg,invtau_qv_nudg,invtau_u_nudg,invtau_v_nudg,             &
     325       du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq,                                    &
     326       dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,tskin,ps,ustar,tkes,uw,vw,q1,q2,       &
     327       orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,          &
     328       heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas)
     329
     330    !program reading forcing of the case study
     331    use netcdf, only: nf90_get_var
     332    implicit none
    335333#include "netcdf.inc"
    336334#include "compar1d.h"
    337335
    338       integer ntime,nlevel,k,t
    339 
    340       real ap(nlevel+1),bp(nlevel+1)
    341       real zz(nlevel,ntime),zzh(nlevel+1)
    342       real pp(nlevel,ntime),pph(nlevel+1)
    343 !profils initiaux
    344       real temp0(nlevel),qv0(nlevel),ql0(nlevel),qi0(nlevel),u0(nlevel),v0(nlevel),tke0(nlevel)
    345       real pp0(nlevel)   
    346       real temp(nlevel,ntime),qv(nlevel,ntime),ql(nlevel,ntime),qi(nlevel,ntime),rh(nlevel,ntime)
    347       real theta(nlevel,ntime),thv(nlevel,ntime),thl(nlevel,ntime),rv(nlevel,ntime)
    348       real u(nlevel,ntime),v(nlevel,ntime),tkes(ntime)
    349       real temp_nudg(nlevel,ntime),qv_nudg(nlevel,ntime),u_nudg(nlevel,ntime),v_nudg(nlevel,ntime)
    350       real invtau_temp_nudg(nlevel,ntime),invtau_qv_nudg(nlevel,ntime),invtau_u_nudg(nlevel,ntime),invtau_v_nudg(nlevel,ntime)
    351       real ug(nlevel,ntime),vg(nlevel,ntime)
    352       real vitw(nlevel,ntime),omega(nlevel,ntime),tke(nlevel,ntime)
    353       real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
    354       real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
    355       real dt(nlevel,ntime),ht(nlevel,ntime),vt(nlevel,ntime)
    356       real dtrad(nlevel,ntime)
    357       real dq(nlevel,ntime),hq(nlevel,ntime),vq(nlevel,ntime)
    358       real dth(nlevel,ntime),hth(nlevel,ntime),vth(nlevel,ntime),hthl(nlevel,ntime)
    359       real dr(nlevel,ntime),hr(nlevel,ntime),vr(nlevel,ntime)
    360       real flat(ntime),sens(ntime),ustar(ntime)
    361       real uw(nlevel,ntime),vw(nlevel,ntime),q1(nlevel,ntime),q2(nlevel,ntime)
    362       real ts(ntime),ps(ntime)
    363       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
    364       real apbp(nlevel+1),resul(nlevel,ntime),resul1(nlevel),resul2(ntime),resul3
    365 
    366 
    367       integer nid, ierr,ierr1,ierr2,rid,i,int_test
    368       integer nbvar3d
    369       parameter(nbvar3d=78)
    370       integer var3didin(nbvar3d),missing_var(nbvar3d)
    371       character*13 name_var(1:nbvar3d)
    372 
    373 
    374 !      data name_var/ &
    375 !     ! coordonnees pression (n+1 niveaux) #4
    376 !     & 'coor_par_a','coor_par_b','height_h','pressure_h',& ! #1-#4
    377 !     ! coordonnees pression (n niveaux) #8
    378 !     &'temp','qv','ql','qi','u','v','tke','pressure',& ! #5-#12
    379 !     ! coordonnees pression + temps #42
    380 !     &'w','omega','ug','vg','uadv','uadvh','uadvv','vadv','vadvh','vadvv','temp_adv','tadvh','tadvv',& !  #13 - #25
    381 !     &'qv_adv','qvadvh','qvadvv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #32
    382 !     & 'radv','radvh','radvv','radcool','q1','q2','ustress','vstress',                           & ! #33 - #40
    383 !     & 'rh','temp_nudging','qv_nudging','u_nudging','v_nudging',                                       & ! #41-45
    384 !     &'height_f','pressure_forc','tempt','theta','thv','thl','qvt','qlt','qit','rv','ut','vt',   & ! #46-58
    385 !     ! coordonnees temps #12
    386 !     &'tkes','sfc_sens_flx','sfc_lat_flx','ts','ps','ustar',&
    387 !     &'orog','albedo','emiss','t_skin','q_skin','mom_rough','heat_rough',&
    388 !     ! scalaires #4
    389 !     &'o3','rugos','clay','sand'/
    390 
    391 
    392 
    393       data name_var/ &
    394      ! coordonnees pression (n+1 niveaux) #4
    395      & 'coor_par_a','coor_par_b','zf','pressure_h',& ! #1-#4
    396      ! coordonnees pression (n niveaux) #8
    397      & 'ta','qv','ql','qi','ua','va','tke','pa',& ! #5-#12
    398      ! coordonnees pression + temps #46
    399      & 'wa','wap','ug','vg','tnua_adv','tnua_advh','tnua_advv','tnva_adv','tnva_advh','tnva_advv','tnta_adv','tnta_advh','tnta_advv',& !  #13 - #25
    400      & 'tnqv_adv','tnqv_advh','tnqv_advv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #32
    401      & 'radv','radvh','radvv','tnta_rad','q1','q2','ustress','vstress',                           & ! #33 - #40
    402      & 'rh','ta_nud','qv_nud','ua_nud','va_nud',                                       & ! #41-45
    403      & 'zh_forc','pa_forc','tat','thetat','thetavt','thetalt','qvt','qlt','qit','rv','uat','vat',   & ! #46-57
    404      & 'nudging_constant_ta', 'nudging_constant_qv', 'nudging_constant_ua', 'nudging_constant_va',           & ! # 58-61
    405      ! coordonnees temps #12
    406      & 'tkes','hfss','hfls','ts_forc','ps_forc','ustar', &                                        ! 62-67
    407      & 'orog','albedo','emiss','t_skin','q_skin','z0','z0h',       &                    ! 68-74
    408      ! scalaires #4                               
    409      & 'O3','rugos','clay','sand'/                                                      ! 75-78
    410 
    411 
    412 !-----------------------------------------------------------------------
    413 ! First check that we are using a version > v2 of the 1D standard format
    414 ! use the difference between 'temp' (old version) and 'ta' (new version)
    415 !-----------------------------------------------------------------------
    416 
    417 
    418          ierr=NF_INQ_VARID(nid,'ta',int_test)
    419          if(ierr/=NF_NOERR) then
    420            print*, '++++++++++++++++++++++++++++++'
    421            print*, 'variable ta missing in cas.nc '
    422            print*, 'You are probably using an obsolete version of the 1D cases'
    423            print*, 'please dowload the last version of the 1D archive from https://lmdz.lmd.jussieu.fr/pub/'
    424            print*, '++++++++++++++++++++++++++++++'
    425            CALL abort_gcm ('mod_1D_cases_read_std','bad version of 1D directory',0)
    426          endif
    427      
    428 !-----------------------------------------------------------------------
    429 ! Checking availability of variable #i in the cas.nc file
    430 !     missing_var=1 if the variable is missing
    431 !-----------------------------------------------------------------------
    432 
    433        do i=1,nbvar3d
    434          missing_var(i)=0.
    435          ierr=NF_INQ_VARID(nid,name_var(i),var3didin(i))
    436          print*, 'name_var(i)', name_var(i), var3didin(i)
    437          if(ierr/=NF_NOERR) then
    438            print *,'Variable manquante dans cas.nc:',i,name_var(i)
    439            ierr=NF_NOERR
    440            missing_var(i)=1
    441          else
    442 
    443 !-----------------------------------------------------------------------
    444 ! Activating keys depending on the presence of specific variables in cas.nc
    445 !-----------------------------------------------------------------------
    446 if ( 1 == 1 ) THEN
    447 ! A MODIFIER: il faudrait dire nudging_temp mais faut le declarer dans compar1d.h etc...       
    448 !           if ( name_var(i) == 'temp_nudging' .and. nint(nudging_t)==0) stop 'Nudging inconsistency temp'
    449             if ( name_var(i) == 'qv_nud' .and. nint(nudging_qv)==0) stop 'Nudging inconsistency qv'
    450             if ( name_var(i) == 'ua_nud' .and. nint(nudging_u)==0) stop 'Nudging inconsistency u'
    451             if ( name_var(i) == 'va_nud' .and. nint(nudging_v)==0) stop 'Nudging inconsistency v'
    452     ELSE
     336    integer ntime,nlevel,k,t
     337
     338    real ap(nlevel+1),bp(nlevel+1)
     339    real zz(nlevel,ntime),zzh(nlevel+1)
     340    real pp(nlevel,ntime),pph(nlevel+1)
     341    !profils initiaux
     342    real temp0(nlevel),qv0(nlevel),ql0(nlevel),qi0(nlevel),u0(nlevel),v0(nlevel),tke0(nlevel)
     343    real pp0(nlevel)   
     344    real temp(nlevel,ntime),qv(nlevel,ntime),ql(nlevel,ntime),qi(nlevel,ntime),rh(nlevel,ntime)
     345    real theta(nlevel,ntime),thv(nlevel,ntime),thl(nlevel,ntime),rv(nlevel,ntime)
     346    real u(nlevel,ntime),v(nlevel,ntime),tkes(ntime)
     347    real temp_nudg(nlevel,ntime),qv_nudg(nlevel,ntime),u_nudg(nlevel,ntime),v_nudg(nlevel,ntime)
     348    real invtau_temp_nudg(nlevel,ntime),invtau_qv_nudg(nlevel,ntime),invtau_u_nudg(nlevel,ntime),invtau_v_nudg(nlevel,ntime)
     349    real ug(nlevel,ntime),vg(nlevel,ntime)
     350    real vitw(nlevel,ntime),omega(nlevel,ntime),tke(nlevel,ntime)
     351    real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
     352    real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
     353    real dt(nlevel,ntime),ht(nlevel,ntime),vt(nlevel,ntime)
     354    real dtrad(nlevel,ntime)
     355    real dq(nlevel,ntime),hq(nlevel,ntime),vq(nlevel,ntime)
     356    real dth(nlevel,ntime),hth(nlevel,ntime),vth(nlevel,ntime),hthl(nlevel,ntime)
     357    real dr(nlevel,ntime),hr(nlevel,ntime),vr(nlevel,ntime)
     358    real flat(ntime),sens(ntime),ustar(ntime)
     359    real uw(nlevel,ntime),vw(nlevel,ntime),q1(nlevel,ntime),q2(nlevel,ntime)
     360    real ts(ntime),tskin(ntime),ps(ntime)
     361    real orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas
     362    real apbp(nlevel+1),resul(nlevel,ntime),resul1(nlevel),resul2(ntime),resul3
     363
     364
     365    integer nid, ierr,ierr1,ierr2,rid,i,int_test
     366    integer nbvar3d
     367    parameter(nbvar3d=78)
     368    integer var3didin(nbvar3d),missing_var(nbvar3d)
     369    character*13 name_var(1:nbvar3d)
     370
     371
     372    !      data name_var/ &
     373    !     ! coordonnees pression (n+1 niveaux) #4
     374    !     & 'coor_par_a','coor_par_b','height_h','pressure_h',& ! #1-#4
     375    !     ! coordonnees pression (n niveaux) #8
     376    !     &'temp','qv','ql','qi','u','v','tke','pressure',& ! #5-#12
     377    !     ! coordonnees pression + temps #42
     378    !     &'w','omega','ug','vg','uadv','uadvh','uadvv','vadv','vadvh','vadvv','temp_adv','tadvh','tadvv',& !  #13 - #25
     379    !     &'qv_adv','qvadvh','qvadvv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #32
     380    !     & 'radv','radvh','radvv','radcool','q1','q2','ustress','vstress',                           & ! #33 - #40
     381    !     & 'rh','temp_nudging','qv_nudging','u_nudging','v_nudging',                                       & ! #41-45
     382    !     &'height_f','pressure_forc','tempt','theta','thv','thl','qvt','qlt','qit','rv','ut','vt',   & ! #46-58
     383    !     ! coordonnees temps #12
     384    !     &'tkes','sfc_sens_flx','sfc_lat_flx','ts','ps','ustar',&
     385    !     &'orog','albedo','emiss','t_skin','q_skin','mom_rough','heat_rough',&
     386    !     ! scalaires #4
     387    !     &'o3','rugos','clay','sand'/
     388
     389
     390
     391    data name_var/ &
     392                                ! coordonnees pression (n+1 niveaux) #4
     393         'coor_par_a','coor_par_b','zf','pressure_h',& ! #1-#4
     394                                ! coordonnees pression (n niveaux) #8
     395         'ta','qv','ql','qi','ua','va','tke','pa',& ! #5-#12
     396                                ! coordonnees pression + temps #46
     397         'wa','wap','ug','vg','tnua_adv','tnua_advh','tnua_advv','tnva_adv','tnva_advh','tnva_advv','tnta_adv','tnta_advh',& !  #13 - #25
     398         'tnta_advv','tnqv_adv','tnqv_advh','tnqv_advv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #32
     399         'radv','radvh','radvv','tnta_rad','q1','q2','ustress','vstress',                           & ! #33 - #40
     400         'rh','ta_nud','qv_nud','ua_nud','va_nud',                                       & ! #41-45
     401         'zh_forc','pa_forc','tat','thetat','thetavt','thetalt','qvt','qlt','qit','rvt','uat','vat',   & ! #46-57
     402         'nudging_constant_ta', 'nudging_constant_qv', 'nudging_constant_ua', 'nudging_constant_va',           & ! # 58-61
     403                                ! coordonnees temps #12
     404         'tkes','hfss','hfls','ts_forc','tskin','ps_forc','ustar', &                     ! 62-68
     405                                  ! scalaires
     406         'orog','albedo','emiss','q_skin','z0','z0h',       &                    ! 69-74
     407         'O3','rugos','clay','sand'/                                                      ! 75-78
     408
     409
     410    !-----------------------------------------------------------------------
     411    ! First check that we are using a version > v2 of the 1D standard format
     412    ! use the difference between 'temp' (old version) and 'ta' (new version)
     413    !-----------------------------------------------------------------------
     414
     415
     416    ierr=NF_INQ_VARID(nid,'ta',int_test)
     417    if(ierr/=NF_NOERR) then
     418       print*, '++++++++++++++++++++++++++++++'
     419       print*, 'variable ta missing in cas.nc '
     420       print*, 'You are probably using an obsolete version of the 1D cases'
     421       print*, 'please dowload the last version of the 1D archive from https://lmdz.lmd.jussieu.fr/pub/'
     422       print*, '++++++++++++++++++++++++++++++'
     423       CALL abort_gcm ('mod_1D_cases_read_std','bad version of 1D directory',0)
     424    endif
     425
     426    !-----------------------------------------------------------------------
     427    ! Checking availability of variable #i in the cas.nc file
     428    !     missing_var=1 if the variable is missing
     429    !-----------------------------------------------------------------------
     430
     431    do i=1,nbvar3d
     432       missing_var(i)=0.
     433       ierr=NF_INQ_VARID(nid,name_var(i),var3didin(i))
     434       print*, 'name_var(i)', name_var(i), var3didin(i)
     435       if(ierr/=NF_NOERR) then
     436          print *,'Variable manquante dans cas.nc:',i,name_var(i)
     437          ierr=NF_NOERR
     438          missing_var(i)=1
     439       else
     440
     441          !-----------------------------------------------------------------------
     442          ! Activating keys depending on the presence of specific variables in cas.nc
     443          !-----------------------------------------------------------------------
     444          if ( 1 == 1 ) THEN
     445             ! A MODIFIER: il faudrait dire nudging_temp mais faut le declarer dans compar1d.h etc...       
     446             !           if ( name_var(i) == 'temp_nudging' .and. nint(nudging_t)==0) stop 'Nudging inconsistency temp'
     447             if ( name_var(i) == 'qv_nud' .and. nint(nudging_qv)==0) stop 'Nudging inconsistency qv'
     448             if ( name_var(i) == 'ua_nud' .and. nint(nudging_u)==0) stop 'Nudging inconsistency u'
     449             if ( name_var(i) == 'va_nud' .and. nint(nudging_v)==0) stop 'Nudging inconsistency v'
     450          ELSE
    453451             print*,'GUIDAGE : CONSISTENCY CHECK DEACTIVATED FOR TESTS of SANDU/REF'
    454     ENDIF
    455 
    456 !-----------------------------------------------------------------------
    457 ! Reading variables 1D (N+1) vertical variables (nlevelp1,lat,lon)
    458 !-----------------------------------------------------------------------
    459            if(i.LE.4) then
    460 #ifdef NC_DOUBLE
    461            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),apbp)
    462 #else
    463            ierr = NF_GET_VAR_REAL(nid,var3didin(i),apbp)
    464 #endif
    465            print *,'read2_cas(apbp), on a lu ',i,name_var(i)
    466            if(ierr/=NF_NOERR) then
    467               print *,'Pb a la lecture de cas.nc: ',name_var(i)
    468               stop "getvarup"
    469            endif
    470 
    471 !-----------------------------------------------------------------------
    472 !  Reading 1D (N) vertical varialbes    (nlevel,lat,lon)   
    473 !-----------------------------------------------------------------------
    474            else if(i.gt.4.and.i.LE.12) then 
    475 #ifdef NC_DOUBLE
    476            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul1)
    477 #else
    478            ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul1)
    479 #endif
    480            print *,'read2_cas(resul1), on a lu ',i,name_var(i)
    481            if(ierr/=NF_NOERR) then
    482               print *,'Pb a la lecture de cas.nc: ',name_var(i)
    483               stop "getvarup"
    484            endif
    485          print*,'Lecture de la variable #i ',i,name_var(i),minval(resul1),maxval(resul1)
    486 
    487 !-----------------------------------------------------------------------
    488 !  Reading 2D tim-vertical variables  (time,nlevel,lat,lon)
    489 !  TBD : seems to be the same as above.
    490 !-----------------------------------------------------------------------
    491            else if(i.gt.12.and.i.LE.61) then
    492 #ifdef NC_DOUBLE
    493            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul)
    494 #else
    495            ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul)
    496 #endif
    497            print *,'read2_cas(resul), on a lu ',i,name_var(i)
    498            if(ierr/=NF_NOERR) then
    499               print *,'Pb a la lecture de cas.nc: ',name_var(i)
    500               stop "getvarup"
    501            endif
    502          print*,'Lecture de la variable #i ',i,name_var(i),minval(resul),maxval(resul)
    503 
    504 !-----------------------------------------------------------------------
    505 !  Reading 1D time variables (time,lat,lon)
    506 !-----------------------------------------------------------------------
    507            else if (i.gt.62.and.i.LE.75) then
    508 #ifdef NC_DOUBLE
    509            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul2)
    510 #else
    511            ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul2)
    512 #endif
    513            print *,'read2_cas(resul2), on a lu ',i,name_var(i)
    514            if(ierr/=NF_NOERR) then
    515               print *,'Pb a la lecture de cas.nc: ',name_var(i)
    516               stop "getvarup"
    517            endif
    518          print*,'Lecture de la variable #i  ',i,name_var(i),minval(resul2),maxval(resul2)
    519 
    520 !-----------------------------------------------------------------------
    521 ! Reading scalar variables (lat,lon)
    522 !-----------------------------------------------------------------------
    523            else
    524 #ifdef NC_DOUBLE
    525            ierr = NF_GET_VAR_DOUBLE(nid,var3didin(i),resul3)
    526 #else
    527            ierr = NF_GET_VAR_REAL(nid,var3didin(i),resul3)
    528 #endif
    529            print *,'read2_cas(resul3), on a lu ',i,name_var(i)
    530            if(ierr/=NF_NOERR) then
    531               print *,'Pb a la lecture de cas.nc: ',name_var(i)
    532               stop "getvarup"
    533            endif
    534          print*,'Lecture de la variable #i ',i,name_var(i),resul3
    535            endif
    536          endif
    537 
    538 !-----------------------------------------------------------------------
    539 ! Attributing variables
    540 !-----------------------------------------------------------------------
    541          select case(i)
    542          !case(1) ; ap=apbp       ! donnees indexees en nlevel+1
    543          ! case(2) ; bp=apbp
    544            case(3) ; zzh=apbp
    545            case(4) ; pph=apbp
    546            case(5) ; temp0=resul1    ! donnees initiales
    547            case(6) ; qv0=resul1
    548            case(7) ; ql0=resul1
    549            case(8) ; qi0=resul1
    550            case(9) ; u0=resul1
    551            case(10) ; v0=resul1
    552            case(11) ; tke0=resul1
    553            case(12) ; pp0=resul1
    554            case(13) ; vitw=resul    ! donnees indexees en nlevel,time
    555            case(14) ; omega=resul
    556            case(15) ; ug=resul
    557            case(16) ; vg=resul
    558            case(17) ; du=resul
    559            case(18) ; hu=resul
    560            case(19) ; vu=resul
    561            case(20) ; dv=resul
    562            case(21) ; hv=resul
    563            case(22) ; vv=resul
    564            case(23) ; dt=resul
    565            case(24) ; ht=resul
    566            case(25) ; vt=resul
    567            case(26) ; dq=resul
    568            case(27) ; hq=resul
    569            case(28) ; vq=resul
    570            case(29) ; dth=resul
    571            case(30) ; hth=resul
    572            case(31) ; vth=resul
    573            case(32) ; hthl=resul
    574            case(33) ; dr=resul
    575            case(34) ; hr=resul
    576            case(35) ; vr=resul
    577            case(36) ; dtrad=resul
    578            case(37) ; q1=resul
    579            case(38) ; q2=resul
    580            case(39) ; uw=resul
    581            case(40) ; vw=resul
    582            case(41) ; rh=resul
    583            case(42) ; temp_nudg=resul
    584            case(43) ; qv_nudg=resul
    585            case(44) ; u_nudg=resul
    586            case(45) ; v_nudg=resul
    587            case(46) ; zz=resul      ! donnees en time,nlevel pour profil initial
    588            case(47) ; pp=resul
    589            case(48) ; temp=resul
    590            case(49) ; theta=resul
    591            case(50) ; thv=resul
    592            case(51) ; thl=resul
    593            case(52) ; qv=resul
    594            case(53) ; ql=resul
    595            case(54) ; qi=resul
    596            case(55) ; rv=resul
    597            case(56) ; u=resul
    598            case(57) ; v=resul
    599            case(58) ; invtau_temp_nudg=resul
    600            case(59) ; invtau_qv_nudg=resul
    601            case(60) ; invtau_u_nudg=resul
    602            case(61) ; invtau_v_nudg=resul
    603            case(62) ; tkes=resul2   ! donnees indexees en time
    604            case(63) ; sens=resul2
    605            case(64) ; flat=resul2
    606            case(65) ; ts=resul2
    607            case(66) ; ps=resul2
    608            case(67) ; ustar=resul2
    609            case(68) ; orog_cas=resul3      ! constantes
    610            case(69) ; albedo_cas=resul3
    611            case(70) ; emiss_cas=resul3
    612            case(71) ; t_skin_cas=resul3
    613            case(72) ; q_skin_cas=resul3
    614            case(73) ; mom_rough=resul3
    615            case(74) ; heat_rough=resul3
    616            case(75) ; o3_cas=resul3       
    617            case(76) ; rugos_cas=resul3
    618            case(77) ; clay_cas=resul3
    619            case(78) ; sand_cas=resul3
    620          end select
    621          resul=0.
    622          resul1=0.
    623          resul2=0.
    624          resul3=0.
     452          ENDIF
     453
     454          !-----------------------------------------------------------------------
     455          ! Reading variables 1D (N+1) vertical variables (nlevelp1,lat,lon)
     456          !-----------------------------------------------------------------------
     457          if(i.LE.4) then
     458             ierr = NF90_GET_VAR(nid,var3didin(i),apbp)
     459             print *,'read_SCM(apbp), on a lu ',i,name_var(i)
     460             if(ierr/=NF_NOERR) then
     461                print *,'B Pb a la lecture de cas.nc: ',name_var(i)
     462                stop "getvarup"
     463             endif
     464
     465             !-----------------------------------------------------------------------
     466             !  Reading 1D (N) vertical varialbes    (nlevel,lat,lon)   
     467             !-----------------------------------------------------------------------
     468          else if(i.gt.4.and.i.LE.12) then 
     469             ierr = NF90_GET_VAR(nid,var3didin(i),resul1)
     470             print *,'read_SCM(resul1), on a lu ',i,name_var(i)
     471             if(ierr/=NF_NOERR) then
     472                print *,'C Pb a la lecture de cas.nc: ',name_var(i)
     473                stop "getvarup"
     474             endif
     475             print*,'Lecture de la variable #i ',i,name_var(i),minval(resul1),maxval(resul1)
     476
     477             !-----------------------------------------------------------------------
     478             !  Reading 2D tim-vertical variables  (time,nlevel,lat,lon)
     479             !  TBD : seems to be the same as above.
     480             !-----------------------------------------------------------------------
     481          else if(i.gt.12.and.i.LE.61) then
     482             ierr = NF90_GET_VAR(nid,var3didin(i),resul)
     483             print *,'read_SCM(resul), on a lu ',i,name_var(i)
     484             if(ierr/=NF_NOERR) then
     485                print *,'D Pb a la lecture de cas.nc: ',name_var(i)
     486                stop "getvarup"
     487             endif
     488             print*,'Lecture de la variable #i ',i,name_var(i),minval(resul),maxval(resul)
     489
     490             !-----------------------------------------------------------------------
     491             !  Reading 1D time variables (time,lat,lon)
     492             !-----------------------------------------------------------------------
     493          else if (i.gt.62.and.i.LE.75) then
     494             ierr = NF90_GET_VAR(nid,var3didin(i),resul2)
     495             print *,'read_SCM(resul2), on a lu ',i,name_var(i)
     496             if(ierr/=NF_NOERR) then
     497                print *,'E Pb a la lecture de cas.nc: ',name_var(i)
     498                stop "getvarup"
     499             endif
     500             print*,'Lecture de la variable #i  ',i,name_var(i),minval(resul2),maxval(resul2)
     501
     502             !-----------------------------------------------------------------------
     503             ! Reading scalar variables (lat,lon)
     504             !-----------------------------------------------------------------------
     505          else
     506             ierr = NF90_GET_VAR(nid,var3didin(i),resul3)
     507             print *,'read_SCM(resul3), on a lu ',i,name_var(i)
     508             if(ierr/=NF_NOERR) then
     509                print *,'F Pb a la lecture de cas.nc: ',name_var(i)
     510                stop "getvarup"
     511             endif
     512             print*,'Lecture de la variable #i ',i,name_var(i),resul3
     513          endif
     514       endif
     515
     516       !-----------------------------------------------------------------------
     517       ! Attributing variables
     518       !-----------------------------------------------------------------------
     519       select case(i)
     520          !case(1) ; ap=apbp       ! donnees indexees en nlevel+1
     521          ! case(2) ; bp=apbp
     522       case(3) ; zzh=apbp
     523       case(4) ; pph=apbp
     524       case(5) ; temp0=resul1    ! donnees initiales
     525       case(6) ; qv0=resul1
     526       case(7) ; ql0=resul1
     527       case(8) ; qi0=resul1
     528       case(9) ; u0=resul1
     529       case(10) ; v0=resul1
     530       case(11) ; tke0=resul1
     531       case(12) ; pp0=resul1
     532       case(13) ; vitw=resul    ! donnees indexees en nlevel,time
     533       case(14) ; omega=resul
     534       case(15) ; ug=resul
     535       case(16) ; vg=resul
     536       case(17) ; du=resul
     537       case(18) ; hu=resul
     538       case(19) ; vu=resul
     539       case(20) ; dv=resul
     540       case(21) ; hv=resul
     541       case(22) ; vv=resul
     542       case(23) ; dt=resul
     543       case(24) ; ht=resul
     544       case(25) ; vt=resul
     545       case(26) ; dq=resul
     546       case(27) ; hq=resul
     547       case(28) ; vq=resul
     548       case(29) ; dth=resul
     549       case(30) ; hth=resul
     550       case(31) ; vth=resul
     551       case(32) ; hthl=resul
     552       case(33) ; dr=resul
     553       case(34) ; hr=resul
     554       case(35) ; vr=resul
     555       case(36) ; dtrad=resul
     556       case(37) ; q1=resul
     557       case(38) ; q2=resul
     558       case(39) ; uw=resul
     559       case(40) ; vw=resul
     560       case(41) ; rh=resul
     561       case(42) ; temp_nudg=resul
     562       case(43) ; qv_nudg=resul
     563       case(44) ; u_nudg=resul
     564       case(45) ; v_nudg=resul
     565       case(46) ; zz=resul      ! donnees en time,nlevel pour profil initial
     566       case(47) ; pp=resul
     567       case(48) ; temp=resul
     568       case(49) ; theta=resul
     569       case(50) ; thv=resul
     570       case(51) ; thl=resul
     571       case(52) ; qv=resul
     572       case(53) ; ql=resul
     573       case(54) ; qi=resul
     574       case(55) ; rv=resul
     575       case(56) ; u=resul
     576       case(57) ; v=resul
     577       case(58) ; invtau_temp_nudg=resul
     578       case(59) ; invtau_qv_nudg=resul
     579       case(60) ; invtau_u_nudg=resul
     580       case(61) ; invtau_v_nudg=resul
     581       case(62) ; tkes=resul2   ! donnees indexees en time
     582       case(63) ; sens=resul2
     583       case(64) ; flat=resul2
     584       case(65) ; ts=resul2
     585       case(66) ; tskin=resul2       
     586       case(67) ; ps=resul2
     587       case(68) ; ustar=resul2
     588       case(69) ; orog_cas=resul3      ! constantes
     589       case(70) ; albedo_cas=resul3
     590       case(71) ; emiss_cas=resul3
     591       case(72) ; q_skin_cas=resul3
     592       case(73) ; mom_rough=resul3
     593       case(74) ; heat_rough=resul3
     594       case(75) ; o3_cas=resul3       
     595       case(76) ; rugos_cas=resul3
     596       case(77) ; clay_cas=resul3
     597       case(78) ; sand_cas=resul3
     598       end select
     599       resul=0.
     600       resul1=0.
     601       resul2=0.
     602       resul3=0.
     603    enddo
     604    print*,'Lecture de la variable APRES ,sens ',minval(sens),maxval(sens)
     605    print*,'Lecture de la variable APRES ,flat ',minval(flat),maxval(flat)
     606
     607    !CR:ATTENTION EN ATTENTE DE REGLER LA QUESTION DU PAS DE TEMPS INITIAL
     608    do t=1,ntime
     609       do k=1,nlevel
     610          temp(k,t)=temp0(k)
     611          qv(k,t)=qv0(k)
     612          ql(k,t)=ql0(k)
     613          qi(k,t)=qi0(k)
     614          u(k,t)=u0(k)
     615          v(k,t)=v0(k)
     616          tke(k,t)=tke0(k)
    625617       enddo
    626          print*,'Lecture de la variable APRES ,sens ',minval(sens),maxval(sens)
    627          print*,'Lecture de la variable APRES ,flat ',minval(flat),maxval(flat)
    628 
    629 !CR:ATTENTION EN ATTENTE DE REGLER LA QUESTION DU PAS DE TEMPS INITIAL
    630        do t=1,ntime
    631           do k=1,nlevel
    632              temp(k,t)=temp0(k)
    633              qv(k,t)=qv0(k)
    634              ql(k,t)=ql0(k)
    635              qi(k,t)=qi0(k)
    636              u(k,t)=u0(k)
    637              v(k,t)=v0(k)
    638              tke(k,t)=tke0(k)
    639           enddo
    640        enddo
    641        !!!! TRAVAIL : EN FONCTION DES DECISIONS SUR LA SPECIFICATION DE W
    642        !!!omega=-vitw*pres*rg/(rd*temp)
    643 !-----------------------------------------------------------------------
    644 
    645          return
    646          END SUBROUTINE read_SCM
    647 !======================================================================
    648 
    649 !======================================================================
    650 
    651 !**********************************************************************************************
    652 
    653 !**********************************************************************************************
    654         SUBROUTINE interp_case_time_std(day,day1,annee_ref                           &
    655 !    &         ,year_cas,day_cas,nt_cas,pdt_forc,nlev_cas                         &
    656      &         ,nt_cas,nlev_cas                                                   &
    657      &         ,ts_cas,ps_cas,plev_cas,t_cas,theta_cas,thv_cas,thl_cas            &
    658      &         ,qv_cas,ql_cas,qi_cas,u_cas,v_cas                                  &
    659      &         ,ug_cas,vg_cas,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas     &
    660      &         ,invtau_temp_nudg_cas,invtau_qv_nudg_cas,invtau_u_nudg_cas,invtau_v_nudg_cas     &
    661      &         ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas             &
    662      &         ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas               &
    663      &         ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas                      &
    664      &         ,lat_cas,sens_cas,ustar_cas                                        &
    665      &         ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                               &
    666 !
    667      &         ,ts_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas               &
    668      &         ,thv_prof_cas,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas     &
    669      &         ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                     &
    670      &         ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas     &
    671      &         ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas     &     
    672      &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas  &
    673      &         ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas                   &
    674      &         ,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas                &
    675      &         ,hq_prof_cas,vq_prof_cas,dth_prof_cas,hth_prof_cas,vth_prof_cas    &
    676      &         ,lat_prof_cas,sens_prof_cas                                        &
    677      &         ,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
    678 
    679 
    680 
    681 
    682 
    683 
    684        implicit none
    685 
    686 !---------------------------------------------------------------------------------------
    687 ! Time interpolation of a 2D field to the timestep corresponding to day
    688 !
    689 ! day: current julian day (e.g. 717538.2)
    690 ! day1: first day of the simulation
    691 ! nt_cas: total nb of data in the forcing
    692 ! pdt_cas: total time interval (in sec) between 2 forcing data
    693 !---------------------------------------------------------------------------------------
     618    enddo
     619!!!! TRAVAIL : EN FONCTION DES DECISIONS SUR LA SPECIFICATION DE W
     620!!!omega=-vitw*pres*rg/(rd*temp)
     621    !-----------------------------------------------------------------------
     622
     623    return
     624  END SUBROUTINE read_SCM
     625  !======================================================================
     626
     627  !======================================================================
     628
     629  !**********************************************************************************************
     630
     631  !**********************************************************************************************
     632  SUBROUTINE interp_case_time_std(day,day1,annee_ref                           &
     633       !    &         ,year_cas,day_cas,nt_cas,pdt_forc,nlev_cas                         &
     634       ,nt_cas,nlev_cas                                                   &
     635       ,ts_cas,tskin_cas,ps_cas,plev_cas,t_cas,theta_cas,thv_cas,thl_cas            &
     636       ,qv_cas,ql_cas,qi_cas,u_cas,v_cas                                  &
     637       ,ug_cas,vg_cas,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas     &
     638       ,invtau_temp_nudg_cas,invtau_qv_nudg_cas,invtau_u_nudg_cas,invtau_v_nudg_cas     &
     639       ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas             &
     640       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas               &
     641       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas                      &
     642       ,lat_cas,sens_cas,ustar_cas                                        &
     643       ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                               &
     644       !
     645       ,ts_prof_cas,tskin_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas               &
     646       ,thv_prof_cas,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas     &
     647       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                     &
     648       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas     &
     649       ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas     &     
     650       ,vitw_prof_cas,omega_prof_cas,tke_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas  &
     651       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas                   &
     652       ,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas                &
     653       ,hq_prof_cas,vq_prof_cas,dth_prof_cas,hth_prof_cas,vth_prof_cas    &
     654       ,lat_prof_cas,sens_prof_cas                                        &
     655       ,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
     656
     657
     658
     659
     660
     661
     662    implicit none
     663
     664    !---------------------------------------------------------------------------------------
     665    ! Time interpolation of a 2D field to the timestep corresponding to day
     666    !
     667    ! day: current julian day (e.g. 717538.2)
     668    ! day1: first day of the simulation
     669    ! nt_cas: total nb of data in the forcing
     670    ! pdt_cas: total time interval (in sec) between 2 forcing data
     671    !---------------------------------------------------------------------------------------
    694672
    695673#include "compar1d.h"
    696674#include "date_cas.h"
    697675
    698 ! inputs:
    699         integer annee_ref
    700         integer nt_cas,nlev_cas
    701         real day, day1,day_cas
    702         real ts_cas(nt_cas),ps_cas(nt_cas)
    703         real plev_cas(nlev_cas,nt_cas)
    704         real t_cas(nlev_cas,nt_cas),theta_cas(nlev_cas,nt_cas)
    705         real thv_cas(nlev_cas,nt_cas), thl_cas(nlev_cas,nt_cas)
    706         real qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas),qi_cas(nlev_cas,nt_cas)
    707         real u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas)
    708         real ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas)
    709         real temp_nudg_cas(nlev_cas,nt_cas),qv_nudg_cas(nlev_cas,nt_cas)
    710         real u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas)
    711 
    712         real invtau_temp_nudg_cas(nlev_cas,nt_cas),invtau_qv_nudg_cas(nlev_cas,nt_cas)
    713         real invtau_u_nudg_cas(nlev_cas,nt_cas),invtau_v_nudg_cas(nlev_cas,nt_cas)
    714 
    715         real vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas),tke_cas(nlev_cas,nt_cas)
    716         real du_cas(nlev_cas,nt_cas),hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas)
    717         real dv_cas(nlev_cas,nt_cas),hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas)
    718         real dt_cas(nlev_cas,nt_cas),ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas)
    719         real dth_cas(nlev_cas,nt_cas),hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas)
    720         real dtrad_cas(nlev_cas,nt_cas)
    721         real dq_cas(nlev_cas,nt_cas),hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas)
    722         real lat_cas(nt_cas),sens_cas(nt_cas),tkes_cas(nt_cas)
    723         real ustar_cas(nt_cas),uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas)
    724         real q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas)
    725 
    726 ! outputs:
    727         real plev_prof_cas(nlev_cas)
    728         real t_prof_cas(nlev_cas),theta_prof_cas(nlev_cas),thl_prof_cas(nlev_cas),thv_prof_cas(nlev_cas)
    729         real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
    730         real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
    731         real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas)
    732         real temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas)
    733         real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
    734 
    735         real invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas)
    736         real invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas)
    737 
    738         real vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas),tke_prof_cas(nlev_cas)
    739         real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
    740         real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
    741         real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas)
    742         real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
    743         real dtrad_prof_cas(nlev_cas)
    744         real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
    745         real lat_prof_cas,sens_prof_cas,tkes_prof_cas,ts_prof_cas,ps_prof_cas,ustar_prof_cas
    746         real uw_prof_cas(nlev_cas),vw_prof_cas(nlev_cas),q1_prof_cas(nlev_cas),q2_prof_cas(nlev_cas)
    747 ! local:
    748         integer it_cas1, it_cas2,k
    749         real timeit,time_cas1,time_cas2,frac
    750 
    751 
    752         print*,'Check time',day1,day_ju_ini_cas,day_deb+1,pdt_cas
    753 !       do k=1,nlev_cas
    754 !       print*,'debut de interp2_case_time, plev_cas=',k,plev_cas(k,1)
    755 !       enddo
    756 
    757 ! On teste si la date du cas AMMA est correcte.
    758 ! C est pour memoire car en fait les fichiers .def
    759 ! sont censes etre corrects.
    760 ! A supprimer a terme (MPL 20150623)
    761 !     if ((forcing_type.eq.10).and.(1.eq.0)) then
    762 ! Check that initial day of the simulation consistent with AMMA case:
    763 !      if (annee_ref.ne.2006) then
    764 !       print*,'Pour AMMA, annee_ref doit etre 2006'
    765 !       print*,'Changer annee_ref dans run.def'
    766 !       stop
    767 !      endif
    768 !      if (annee_ref.eq.2006 .and. day1.lt.day_cas) then
    769 !       print*,'AMMA a debute le 10 juillet 2006',day1,day_cas
    770 !       print*,'Changer dayref dans run.def'
    771 !       stop
    772 !      endif
    773 !      if (annee_ref.eq.2006 .and. day1.gt.day_cas+1) then
    774 !       print*,'AMMA a fini le 11 juillet'
    775 !       print*,'Changer dayref ou nday dans run.def'
    776 !       stop
    777 !      endif
    778 !      endif
    779 
    780 ! Determine timestep relative to the 1st day:
    781 !       timeit=(day-day1)*86400.
    782 !       if (annee_ref.eq.1992) then
    783 !        timeit=(day-day_cas)*86400.
    784 !       else
    785 !        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
    786 !       endif
    787       timeit=(day-day_ju_ini_cas)*86400
    788       print *,'day=',day
    789       print *,'day_ju_ini_cas=',day_ju_ini_cas
    790       print *,'pdt_cas=',pdt_cas
    791       print *,'timeit=',timeit
    792       print *,'nt_cas=',nt_cas
    793 
    794 ! Determine the closest observation times:
    795 !       it_cas1=INT(timeit/pdt_cas)+1
    796 !       it_cas2=it_cas1 + 1
    797 !       time_cas1=(it_cas1-1)*pdt_cas
    798 !       time_cas2=(it_cas2-1)*pdt_cas
    799 
    800        it_cas1=INT(timeit/pdt_cas)+1
    801        IF (it_cas1 .EQ. nt_cas) THEN
     676    ! inputs:
     677    integer annee_ref
     678    integer nt_cas,nlev_cas
     679    real day, day1,day_cas
     680    real ts_cas(nt_cas),tskin_cas(nt_cas),ps_cas(nt_cas)
     681    real plev_cas(nlev_cas,nt_cas)
     682    real t_cas(nlev_cas,nt_cas),theta_cas(nlev_cas,nt_cas)
     683    real thv_cas(nlev_cas,nt_cas), thl_cas(nlev_cas,nt_cas)
     684    real qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas),qi_cas(nlev_cas,nt_cas)
     685    real u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas)
     686    real ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas)
     687    real temp_nudg_cas(nlev_cas,nt_cas),qv_nudg_cas(nlev_cas,nt_cas)
     688    real u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas)
     689
     690    real invtau_temp_nudg_cas(nlev_cas,nt_cas),invtau_qv_nudg_cas(nlev_cas,nt_cas)
     691    real invtau_u_nudg_cas(nlev_cas,nt_cas),invtau_v_nudg_cas(nlev_cas,nt_cas)
     692
     693    real vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas),tke_cas(nlev_cas,nt_cas)
     694    real du_cas(nlev_cas,nt_cas),hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas)
     695    real dv_cas(nlev_cas,nt_cas),hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas)
     696    real dt_cas(nlev_cas,nt_cas),ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas)
     697    real dth_cas(nlev_cas,nt_cas),hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas)
     698    real dtrad_cas(nlev_cas,nt_cas)
     699    real dq_cas(nlev_cas,nt_cas),hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas)
     700    real lat_cas(nt_cas),sens_cas(nt_cas),tkes_cas(nt_cas)
     701    real ustar_cas(nt_cas),uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas)
     702    real q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas)
     703
     704    ! outputs:
     705    real plev_prof_cas(nlev_cas)
     706    real t_prof_cas(nlev_cas),theta_prof_cas(nlev_cas),thl_prof_cas(nlev_cas),thv_prof_cas(nlev_cas)
     707    real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
     708    real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
     709    real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas)
     710    real temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas)
     711    real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
     712
     713    real invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas)
     714    real invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas)
     715
     716    real vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas),tke_prof_cas(nlev_cas)
     717    real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
     718    real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
     719    real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas)
     720    real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
     721    real dtrad_prof_cas(nlev_cas)
     722    real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
     723    real lat_prof_cas,sens_prof_cas,tkes_prof_cas,ts_prof_cas,tskin_prof_cas,ps_prof_cas,ustar_prof_cas
     724    real uw_prof_cas(nlev_cas),vw_prof_cas(nlev_cas),q1_prof_cas(nlev_cas),q2_prof_cas(nlev_cas)
     725    ! local:
     726    integer it_cas1, it_cas2,k
     727    real timeit,time_cas1,time_cas2,frac
     728
     729
     730    print*,'Check time',day1,day_ju_ini_cas,day_deb+1,pdt_cas
     731    !       do k=1,nlev_cas
     732    !       print*,'debut de interp2_case_time, plev_cas=',k,plev_cas(k,1)
     733    !       enddo
     734
     735    ! On teste si la date du cas AMMA est correcte.
     736    ! C est pour memoire car en fait les fichiers .def
     737    ! sont censes etre corrects.
     738    ! A supprimer a terme (MPL 20150623)
     739    !     if ((forcing_type.eq.10).and.(1.eq.0)) then
     740    ! Check that initial day of the simulation consistent with AMMA case:
     741    !      if (annee_ref.ne.2006) then
     742    !       print*,'Pour AMMA, annee_ref doit etre 2006'
     743    !       print*,'Changer annee_ref dans run.def'
     744    !       stop
     745    !      endif
     746    !      if (annee_ref.eq.2006 .and. day1.lt.day_cas) then
     747    !       print*,'AMMA a debute le 10 juillet 2006',day1,day_cas
     748    !       print*,'Changer dayref dans run.def'
     749    !       stop
     750    !      endif
     751    !      if (annee_ref.eq.2006 .and. day1.gt.day_cas+1) then
     752    !       print*,'AMMA a fini le 11 juillet'
     753    !       print*,'Changer dayref ou nday dans run.def'
     754    !       stop
     755    !      endif
     756    !      endif
     757
     758    ! Determine timestep relative to the 1st day:
     759    !       timeit=(day-day1)*86400.
     760    !       if (annee_ref.eq.1992) then
     761    !        timeit=(day-day_cas)*86400.
     762    !       else
     763    !        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
     764    !       endif
     765    timeit=(day-day_ju_ini_cas)*86400
     766    print *,'day=',day
     767    print *,'day_ju_ini_cas=',day_ju_ini_cas
     768    print *,'pdt_cas=',pdt_cas
     769    print *,'timeit=',timeit
     770    print *,'nt_cas=',nt_cas
     771
     772    ! Determine the closest observation times:
     773    !       it_cas1=INT(timeit/pdt_cas)+1
     774    !       it_cas2=it_cas1 + 1
     775    !       time_cas1=(it_cas1-1)*pdt_cas
     776    !       time_cas2=(it_cas2-1)*pdt_cas
     777
     778    it_cas1=INT(timeit/pdt_cas)+1
     779    IF (it_cas1 .EQ. nt_cas) THEN
    802780       it_cas2=it_cas1
    803        ELSE
     781    ELSE
    804782       it_cas2=it_cas1 + 1
    805        ENDIF
    806        time_cas1=(it_cas1-1)*pdt_cas
    807        time_cas2=(it_cas2-1)*pdt_cas
    808 !     print *,'timeit,pdt_cas,nt_cas=',timeit,pdt_cas,nt_cas
    809 !     print *,'it_cas1,it_cas2,time_cas1,time_cas2=',it_cas1,it_cas2,time_cas1,time_cas2
    810 
    811        if (it_cas1 .gt. nt_cas) then
    812         write(*,*) 'PB-stop: day, day_ju_ini_cas,it_cas1, it_cas2, timeit: '            &
    813      &        ,day,day_ju_ini_cas,it_cas1,it_cas2,timeit
    814         stop
    815        endif
    816 
    817 ! time interpolation:
    818        IF (it_cas1 .EQ. it_cas2) THEN
    819           frac=0.
    820        ELSE
    821           frac=(time_cas2-timeit)/(time_cas2-time_cas1)
    822           frac=max(frac,0.0)
    823        ENDIF
    824 
    825        lat_prof_cas = lat_cas(it_cas2)                                   &
    826      &          -frac*(lat_cas(it_cas2)-lat_cas(it_cas1))
    827        sens_prof_cas = sens_cas(it_cas2)                                 &
    828      &          -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
    829        tkes_prof_cas = tkes_cas(it_cas2)                                   &
    830      &          -frac*(tkes_cas(it_cas2)-tkes_cas(it_cas1))
    831        ts_prof_cas = ts_cas(it_cas2)                                     &
    832      &          -frac*(ts_cas(it_cas2)-ts_cas(it_cas1))
    833        ps_prof_cas = ps_cas(it_cas2)                                     &
    834      &          -frac*(ps_cas(it_cas2)-ps_cas(it_cas1))
    835        ustar_prof_cas = ustar_cas(it_cas2)                               &
    836      &          -frac*(ustar_cas(it_cas2)-ustar_cas(it_cas1))
    837 
    838        do k=1,nlev_cas
    839         plev_prof_cas(k) = plev_cas(k,it_cas2)                           &     
    840      &          -frac*(plev_cas(k,it_cas2)-plev_cas(k,it_cas1))
    841         t_prof_cas(k) = t_cas(k,it_cas2)                                 &       
    842      &          -frac*(t_cas(k,it_cas2)-t_cas(k,it_cas1))
    843         !print *,'k,frac,plev_cas1,plev_cas2=',k,frac,plev_cas(k,it_cas1),plev_cas(k,it_cas2)
    844         theta_prof_cas(k) = theta_cas(k,it_cas2)                         &                     
    845      &          -frac*(theta_cas(k,it_cas2)-theta_cas(k,it_cas1))
    846         thv_prof_cas(k) = thv_cas(k,it_cas2)                             &         
    847      &          -frac*(thv_cas(k,it_cas2)-thv_cas(k,it_cas1))
    848         thl_prof_cas(k) = thl_cas(k,it_cas2)                             &             
    849      &          -frac*(thl_cas(k,it_cas2)-thl_cas(k,it_cas1))
    850         qv_prof_cas(k) = qv_cas(k,it_cas2)                               &
    851      &          -frac*(qv_cas(k,it_cas2)-qv_cas(k,it_cas1))
    852         ql_prof_cas(k) = ql_cas(k,it_cas2)                               &
    853      &          -frac*(ql_cas(k,it_cas2)-ql_cas(k,it_cas1))
    854         qi_prof_cas(k) = qi_cas(k,it_cas2)                               &
    855      &          -frac*(qi_cas(k,it_cas2)-qi_cas(k,it_cas1))
    856         u_prof_cas(k) = u_cas(k,it_cas2)                                 &
    857      &          -frac*(u_cas(k,it_cas2)-u_cas(k,it_cas1))
    858         v_prof_cas(k) = v_cas(k,it_cas2)                                 &
    859      &          -frac*(v_cas(k,it_cas2)-v_cas(k,it_cas1))
    860         ug_prof_cas(k) = ug_cas(k,it_cas2)                               &
    861      &          -frac*(ug_cas(k,it_cas2)-ug_cas(k,it_cas1))
    862         vg_prof_cas(k) = vg_cas(k,it_cas2)                               &
    863      &          -frac*(vg_cas(k,it_cas2)-vg_cas(k,it_cas1))
    864         temp_nudg_prof_cas(k) = temp_nudg_cas(k,it_cas2)                    &
    865      &          -frac*(temp_nudg_cas(k,it_cas2)-temp_nudg_cas(k,it_cas1))
    866         qv_nudg_prof_cas(k) = qv_nudg_cas(k,it_cas2)                        &
    867      &          -frac*(qv_nudg_cas(k,it_cas2)-qv_nudg_cas(k,it_cas1))
    868         u_nudg_prof_cas(k) = u_nudg_cas(k,it_cas2)                          &
    869      &          -frac*(u_nudg_cas(k,it_cas2)-u_nudg_cas(k,it_cas1))
    870         v_nudg_prof_cas(k) = v_nudg_cas(k,it_cas2)                          &
    871      &          -frac*(v_nudg_cas(k,it_cas2)-v_nudg_cas(k,it_cas1))
    872         invtau_temp_nudg_prof_cas(k) = invtau_temp_nudg_cas(k,it_cas2)                    &
    873      &          -frac*(invtau_temp_nudg_cas(k,it_cas2)-invtau_temp_nudg_cas(k,it_cas1))
    874         invtau_qv_nudg_prof_cas(k) = invtau_qv_nudg_cas(k,it_cas2)                        &
    875      &          -frac*(invtau_qv_nudg_cas(k,it_cas2)-invtau_qv_nudg_cas(k,it_cas1))
    876         invtau_u_nudg_prof_cas(k) = invtau_u_nudg_cas(k,it_cas2)                          &
    877      &          -frac*(invtau_u_nudg_cas(k,it_cas2)-invtau_u_nudg_cas(k,it_cas1))
    878         invtau_v_nudg_prof_cas(k) = invtau_v_nudg_cas(k,it_cas2)                          &
    879      &          -frac*(invtau_v_nudg_cas(k,it_cas2)-invtau_v_nudg_cas(k,it_cas1))
    880         vitw_prof_cas(k) = vitw_cas(k,it_cas2)                           &
    881      &          -frac*(vitw_cas(k,it_cas2)-vitw_cas(k,it_cas1))
    882         omega_prof_cas(k) = omega_cas(k,it_cas2)                         &
    883      &          -frac*(omega_cas(k,it_cas2)-omega_cas(k,it_cas1))
    884         tke_prof_cas(k) = tke_cas(k,it_cas2)                         &
    885      &          -frac*(tke_cas(k,it_cas2)-tke_cas(k,it_cas1))
    886         du_prof_cas(k) = du_cas(k,it_cas2)                               &
    887      &          -frac*(du_cas(k,it_cas2)-du_cas(k,it_cas1))
    888         hu_prof_cas(k) = hu_cas(k,it_cas2)                               &
    889      &          -frac*(hu_cas(k,it_cas2)-hu_cas(k,it_cas1))
    890         vu_prof_cas(k) = vu_cas(k,it_cas2)                               &
    891      &          -frac*(vu_cas(k,it_cas2)-vu_cas(k,it_cas1))
    892         dv_prof_cas(k) = dv_cas(k,it_cas2)                               &
    893      &          -frac*(dv_cas(k,it_cas2)-dv_cas(k,it_cas1))
    894         hv_prof_cas(k) = hv_cas(k,it_cas2)                               &
    895      &          -frac*(hv_cas(k,it_cas2)-hv_cas(k,it_cas1))
    896         vv_prof_cas(k) = vv_cas(k,it_cas2)                               &
    897      &          -frac*(vv_cas(k,it_cas2)-vv_cas(k,it_cas1))
    898         dt_prof_cas(k) = dt_cas(k,it_cas2)                               &
    899      &          -frac*(dt_cas(k,it_cas2)-dt_cas(k,it_cas1))
    900         ht_prof_cas(k) = ht_cas(k,it_cas2)                               &
    901      &          -frac*(ht_cas(k,it_cas2)-ht_cas(k,it_cas1))
    902         vt_prof_cas(k) = vt_cas(k,it_cas2)                               &
    903      &          -frac*(vt_cas(k,it_cas2)-vt_cas(k,it_cas1))
    904         dth_prof_cas(k) = dth_cas(k,it_cas2)                             &
    905      &          -frac*(dth_cas(k,it_cas2)-dth_cas(k,it_cas1))
    906         hth_prof_cas(k) = hth_cas(k,it_cas2)                             &
    907      &          -frac*(hth_cas(k,it_cas2)-hth_cas(k,it_cas1))
    908         vth_prof_cas(k) = vth_cas(k,it_cas2)                             &
    909      &          -frac*(vth_cas(k,it_cas2)-vth_cas(k,it_cas1))
    910         dtrad_prof_cas(k) = dtrad_cas(k,it_cas2)                         &
    911      &          -frac*(dtrad_cas(k,it_cas2)-dtrad_cas(k,it_cas1))
    912         dq_prof_cas(k) = dq_cas(k,it_cas2)                               &
    913      &          -frac*(dq_cas(k,it_cas2)-dq_cas(k,it_cas1))
    914         hq_prof_cas(k) = hq_cas(k,it_cas2)                               &
    915      &          -frac*(hq_cas(k,it_cas2)-hq_cas(k,it_cas1))
    916         vq_prof_cas(k) = vq_cas(k,it_cas2)                               &
    917      &          -frac*(vq_cas(k,it_cas2)-vq_cas(k,it_cas1))
     783    ENDIF
     784    time_cas1=(it_cas1-1)*pdt_cas
     785    time_cas2=(it_cas2-1)*pdt_cas
     786    !     print *,'timeit,pdt_cas,nt_cas=',timeit,pdt_cas,nt_cas
     787    !     print *,'it_cas1,it_cas2,time_cas1,time_cas2=',it_cas1,it_cas2,time_cas1,time_cas2
     788
     789    if (it_cas1 .gt. nt_cas) then
     790       write(*,*) 'PB-stop: day, day_ju_ini_cas,it_cas1, it_cas2, timeit: '            &
     791            ,day,day_ju_ini_cas,it_cas1,it_cas2,timeit
     792       stop
     793    endif
     794
     795    ! time interpolation:
     796    IF (it_cas1 .EQ. it_cas2) THEN
     797       frac=0.
     798    ELSE
     799       frac=(time_cas2-timeit)/(time_cas2-time_cas1)
     800       frac=max(frac,0.0)
     801    ENDIF
     802
     803    lat_prof_cas = lat_cas(it_cas2)                                   &
     804         -frac*(lat_cas(it_cas2)-lat_cas(it_cas1))
     805    sens_prof_cas = sens_cas(it_cas2)                                 &
     806         -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
     807    tkes_prof_cas = tkes_cas(it_cas2)                                   &
     808         -frac*(tkes_cas(it_cas2)-tkes_cas(it_cas1))
     809    ts_prof_cas = ts_cas(it_cas2)                                     &
     810         -frac*(ts_cas(it_cas2)-ts_cas(it_cas1))
     811    tskin_prof_cas = tskin_cas(it_cas2)                                     &
     812         -frac*(tskin_cas(it_cas2)-tskin_cas(it_cas1))
     813    ps_prof_cas = ps_cas(it_cas2)                                     &
     814         -frac*(ps_cas(it_cas2)-ps_cas(it_cas1))
     815    ustar_prof_cas = ustar_cas(it_cas2)                               &
     816         -frac*(ustar_cas(it_cas2)-ustar_cas(it_cas1))
     817
     818    do k=1,nlev_cas
     819       plev_prof_cas(k) = plev_cas(k,it_cas2)                           &     
     820            -frac*(plev_cas(k,it_cas2)-plev_cas(k,it_cas1))
     821       t_prof_cas(k) = t_cas(k,it_cas2)                                 &       
     822            -frac*(t_cas(k,it_cas2)-t_cas(k,it_cas1))
     823       !print *,'k,frac,plev_cas1,plev_cas2=',k,frac,plev_cas(k,it_cas1),plev_cas(k,it_cas2)
     824       theta_prof_cas(k) = theta_cas(k,it_cas2)                         &                     
     825            -frac*(theta_cas(k,it_cas2)-theta_cas(k,it_cas1))
     826       thv_prof_cas(k) = thv_cas(k,it_cas2)                             &         
     827            -frac*(thv_cas(k,it_cas2)-thv_cas(k,it_cas1))
     828       thl_prof_cas(k) = thl_cas(k,it_cas2)                             &             
     829            -frac*(thl_cas(k,it_cas2)-thl_cas(k,it_cas1))
     830       qv_prof_cas(k) = qv_cas(k,it_cas2)                               &
     831            -frac*(qv_cas(k,it_cas2)-qv_cas(k,it_cas1))
     832       ql_prof_cas(k) = ql_cas(k,it_cas2)                               &
     833            -frac*(ql_cas(k,it_cas2)-ql_cas(k,it_cas1))
     834       qi_prof_cas(k) = qi_cas(k,it_cas2)                               &
     835            -frac*(qi_cas(k,it_cas2)-qi_cas(k,it_cas1))
     836       u_prof_cas(k) = u_cas(k,it_cas2)                                 &
     837            -frac*(u_cas(k,it_cas2)-u_cas(k,it_cas1))
     838       v_prof_cas(k) = v_cas(k,it_cas2)                                 &
     839            -frac*(v_cas(k,it_cas2)-v_cas(k,it_cas1))
     840       ug_prof_cas(k) = ug_cas(k,it_cas2)                               &
     841            -frac*(ug_cas(k,it_cas2)-ug_cas(k,it_cas1))
     842       vg_prof_cas(k) = vg_cas(k,it_cas2)                               &
     843            -frac*(vg_cas(k,it_cas2)-vg_cas(k,it_cas1))
     844       temp_nudg_prof_cas(k) = temp_nudg_cas(k,it_cas2)                    &
     845            -frac*(temp_nudg_cas(k,it_cas2)-temp_nudg_cas(k,it_cas1))
     846       qv_nudg_prof_cas(k) = qv_nudg_cas(k,it_cas2)                        &
     847            -frac*(qv_nudg_cas(k,it_cas2)-qv_nudg_cas(k,it_cas1))
     848       u_nudg_prof_cas(k) = u_nudg_cas(k,it_cas2)                          &
     849            -frac*(u_nudg_cas(k,it_cas2)-u_nudg_cas(k,it_cas1))
     850       v_nudg_prof_cas(k) = v_nudg_cas(k,it_cas2)                          &
     851            -frac*(v_nudg_cas(k,it_cas2)-v_nudg_cas(k,it_cas1))
     852       invtau_temp_nudg_prof_cas(k) = invtau_temp_nudg_cas(k,it_cas2)                    &
     853            -frac*(invtau_temp_nudg_cas(k,it_cas2)-invtau_temp_nudg_cas(k,it_cas1))
     854       invtau_qv_nudg_prof_cas(k) = invtau_qv_nudg_cas(k,it_cas2)                        &
     855            -frac*(invtau_qv_nudg_cas(k,it_cas2)-invtau_qv_nudg_cas(k,it_cas1))
     856       invtau_u_nudg_prof_cas(k) = invtau_u_nudg_cas(k,it_cas2)                          &
     857            -frac*(invtau_u_nudg_cas(k,it_cas2)-invtau_u_nudg_cas(k,it_cas1))
     858       invtau_v_nudg_prof_cas(k) = invtau_v_nudg_cas(k,it_cas2)                          &
     859            -frac*(invtau_v_nudg_cas(k,it_cas2)-invtau_v_nudg_cas(k,it_cas1))
     860       vitw_prof_cas(k) = vitw_cas(k,it_cas2)                           &
     861            -frac*(vitw_cas(k,it_cas2)-vitw_cas(k,it_cas1))
     862       omega_prof_cas(k) = omega_cas(k,it_cas2)                         &
     863            -frac*(omega_cas(k,it_cas2)-omega_cas(k,it_cas1))
     864       tke_prof_cas(k) = tke_cas(k,it_cas2)                         &
     865            -frac*(tke_cas(k,it_cas2)-tke_cas(k,it_cas1))
     866       du_prof_cas(k) = du_cas(k,it_cas2)                               &
     867            -frac*(du_cas(k,it_cas2)-du_cas(k,it_cas1))
     868       hu_prof_cas(k) = hu_cas(k,it_cas2)                               &
     869            -frac*(hu_cas(k,it_cas2)-hu_cas(k,it_cas1))
     870       vu_prof_cas(k) = vu_cas(k,it_cas2)                               &
     871            -frac*(vu_cas(k,it_cas2)-vu_cas(k,it_cas1))
     872       dv_prof_cas(k) = dv_cas(k,it_cas2)                               &
     873            -frac*(dv_cas(k,it_cas2)-dv_cas(k,it_cas1))
     874       hv_prof_cas(k) = hv_cas(k,it_cas2)                               &
     875            -frac*(hv_cas(k,it_cas2)-hv_cas(k,it_cas1))
     876       vv_prof_cas(k) = vv_cas(k,it_cas2)                               &
     877            -frac*(vv_cas(k,it_cas2)-vv_cas(k,it_cas1))
     878       dt_prof_cas(k) = dt_cas(k,it_cas2)                               &
     879            -frac*(dt_cas(k,it_cas2)-dt_cas(k,it_cas1))
     880       ht_prof_cas(k) = ht_cas(k,it_cas2)                               &
     881            -frac*(ht_cas(k,it_cas2)-ht_cas(k,it_cas1))
     882       vt_prof_cas(k) = vt_cas(k,it_cas2)                               &
     883            -frac*(vt_cas(k,it_cas2)-vt_cas(k,it_cas1))
     884       dth_prof_cas(k) = dth_cas(k,it_cas2)                             &
     885            -frac*(dth_cas(k,it_cas2)-dth_cas(k,it_cas1))
     886       hth_prof_cas(k) = hth_cas(k,it_cas2)                             &
     887            -frac*(hth_cas(k,it_cas2)-hth_cas(k,it_cas1))
     888       vth_prof_cas(k) = vth_cas(k,it_cas2)                             &
     889            -frac*(vth_cas(k,it_cas2)-vth_cas(k,it_cas1))
     890       dtrad_prof_cas(k) = dtrad_cas(k,it_cas2)                         &
     891            -frac*(dtrad_cas(k,it_cas2)-dtrad_cas(k,it_cas1))
     892       dq_prof_cas(k) = dq_cas(k,it_cas2)                               &
     893            -frac*(dq_cas(k,it_cas2)-dq_cas(k,it_cas1))
     894       hq_prof_cas(k) = hq_cas(k,it_cas2)                               &
     895            -frac*(hq_cas(k,it_cas2)-hq_cas(k,it_cas1))
     896       vq_prof_cas(k) = vq_cas(k,it_cas2)                               &
     897            -frac*(vq_cas(k,it_cas2)-vq_cas(k,it_cas1))
    918898       uw_prof_cas(k) = uw_cas(k,it_cas2)                                &
    919      &          -frac*(uw_cas(k,it_cas2)-uw_cas(k,it_cas1))
     899            -frac*(uw_cas(k,it_cas2)-uw_cas(k,it_cas1))
    920900       vw_prof_cas(k) = vw_cas(k,it_cas2)                                &
    921      &          -frac*(vw_cas(k,it_cas2)-vw_cas(k,it_cas1))
     901            -frac*(vw_cas(k,it_cas2)-vw_cas(k,it_cas1))
    922902       q1_prof_cas(k) = q1_cas(k,it_cas2)                                &
    923      &          -frac*(q1_cas(k,it_cas2)-q1_cas(k,it_cas1))
     903            -frac*(q1_cas(k,it_cas2)-q1_cas(k,it_cas1))
    924904       q2_prof_cas(k) = q2_cas(k,it_cas2)                                &
    925      &          -frac*(q2_cas(k,it_cas2)-q2_cas(k,it_cas1))
    926         enddo
    927 
    928         return
    929         END SUBROUTINE interp_case_time_std
    930 
    931 !**********************************************************************************************
    932 !=====================================================================
    933        SUBROUTINE interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                           &
    934      &         ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas                                       &
    935      &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                              &
    936      &         ,ug_prof_cas,vg_prof_cas                                                                &
    937      &         ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                    &
    938      &         ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas &     
    939      &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                              &
    940      &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                &
    941      &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
    942      &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                 &
    943 !
    944      &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas                                        &
    945      &         ,qv_mod_cas,ql_mod_cas,qi_mod_cas,u_mod_cas,v_mod_cas                                   &
    946      &         ,ug_mod_cas,vg_mod_cas                                                                  &
    947      &         ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                        &
    948      &         ,invtau_temp_nudg_mod_cas,invtau_qv_nudg_mod_cas,invtau_u_nudg_mod_cas,invtau_v_nudg_mod_cas                        &     
    949      &         ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                    &
    950      &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                      &
    951      &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas        &
    952      &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
    953  
    954        implicit none
    955  
     905            -frac*(q2_cas(k,it_cas2)-q2_cas(k,it_cas1))
     906    enddo
     907
     908    return
     909  END SUBROUTINE interp_case_time_std
     910
     911  !**********************************************************************************************
     912  !=====================================================================
     913  SUBROUTINE interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                           &
     914       ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas                                       &
     915       ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                              &
     916       ,ug_prof_cas,vg_prof_cas                                                                &
     917       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                    &
     918       ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas &     
     919       ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                              &
     920       ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                &
     921       ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
     922       ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                 &
     923       !
     924       ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas                                        &
     925       ,qv_mod_cas,ql_mod_cas,qi_mod_cas,u_mod_cas,v_mod_cas                                   &
     926       ,ug_mod_cas,vg_mod_cas                                                                  &
     927       ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                        &
     928       ,invtau_temp_nudg_mod_cas,invtau_qv_nudg_mod_cas,invtau_u_nudg_mod_cas,invtau_v_nudg_mod_cas                        &     
     929       ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                    &
     930       ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                      &
     931       ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas        &
     932       ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
     933
     934    implicit none
     935
    956936#include "YOMCST.h"
    957937#include "dimensions.h"
    958938
    959 !-------------------------------------------------------------------------
    960 ! Vertical interpolation of generic case forcing data onto mod_casel levels
    961 !-------------------------------------------------------------------------
    962  
    963        integer nlevmax
    964        parameter (nlevmax=41)
    965        integer nlev_cas,mxcalc
    966 !       real play(llm), plev_prof(nlevmax)
    967 !       real t_prof(nlevmax),q_prof(nlevmax)
    968 !       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
    969 !       real ht_prof(nlevmax),vt_prof(nlevmax)
    970 !       real hq_prof(nlevmax),vq_prof(nlevmax)
    971  
    972        real play(llm), plev(llm+1), plev_prof_cas(nlev_cas)
    973        real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
    974        real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
    975        real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
    976        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)
    977        real temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas)
    978        real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
    979        real invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas)
    980        real invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas)
    981 
    982        real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
    983        real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
    984        real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
    985        real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
    986        real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
    987  
    988        real t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm)
    989        real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
    990        real u_mod_cas(llm),v_mod_cas(llm)
    991        real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm),tke_mod_cas(llm+1)
    992        real temp_nudg_mod_cas(llm),qv_nudg_mod_cas(llm)
    993        real u_nudg_mod_cas(llm),v_nudg_mod_cas(llm)
    994        real invtau_temp_nudg_mod_cas(llm),invtau_qv_nudg_mod_cas(llm)
    995        real invtau_u_nudg_mod_cas(llm),invtau_v_nudg_mod_cas(llm)
    996        real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
    997        real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
    998        real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
    999        real dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm)
    1000        real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
    1001  
    1002        integer l,k,k1,k2
    1003        real frac,frac1,frac2,fact
    1004  
    1005 
    1006 
    1007 ! for variables defined at the middle of layers
    1008 
    1009        do l = 1, llm
    1010 
    1011         if (play(l).ge.plev_prof_cas(nlev_cas)) then
    1012  
    1013         mxcalc=l
    1014 !        print *,'debut interp2, mxcalc=',mxcalc
    1015          k1=0
    1016          k2=0
    1017 
    1018          if (play(l).le.plev_prof_cas(1)) then
    1019 
    1020          do k = 1, nlev_cas-1
    1021           if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
    1022             k1=k
    1023             k2=k+1
    1024           endif
    1025          enddo
    1026 
    1027          if (k1.eq.0 .or. k2.eq.0) then
    1028           write(*,*) 'PB! k1, k2 = ',k1,k2
    1029           write(*,*) 'l,play(l) = ',l,play(l)/100
    1030          do k = 1, nlev_cas-1
    1031           write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
    1032          enddo
    1033          endif
    1034 
    1035 
    1036 
    1037          frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
    1038          
    1039          t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
    1040          theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
    1041          if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
    1042          thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1))
    1043          thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1))
    1044          qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1))
    1045          ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1))
    1046          qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1))
    1047          u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
    1048          v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
    1049          ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
    1050          vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
    1051          temp_nudg_mod_cas(l)= temp_nudg_prof_cas(k2) - frac*(temp_nudg_prof_cas(k2)-temp_nudg_prof_cas(k1))
    1052          qv_nudg_mod_cas(l)= qv_nudg_prof_cas(k2) - frac*(qv_nudg_prof_cas(k2)-qv_nudg_prof_cas(k1))
    1053          u_nudg_mod_cas(l)= u_nudg_prof_cas(k2) - frac*(u_nudg_prof_cas(k2)-u_nudg_prof_cas(k1))
    1054          v_nudg_mod_cas(l)= v_nudg_prof_cas(k2) - frac*(v_nudg_prof_cas(k2)-v_nudg_prof_cas(k1))
    1055 
    1056          invtau_temp_nudg_mod_cas(l)= invtau_temp_nudg_prof_cas(k2) - frac*(invtau_temp_nudg_prof_cas(k2)-invtau_temp_nudg_prof_cas(k1))
    1057          invtau_qv_nudg_mod_cas(l)= invtau_qv_nudg_prof_cas(k2) - frac*(invtau_qv_nudg_prof_cas(k2)-invtau_qv_nudg_prof_cas(k1))
    1058          invtau_u_nudg_mod_cas(l)= invtau_u_nudg_prof_cas(k2) - frac*(invtau_u_nudg_prof_cas(k2)-invtau_u_nudg_prof_cas(k1))
    1059          invtau_v_nudg_mod_cas(l)= invtau_v_nudg_prof_cas(k2) - frac*(invtau_v_nudg_prof_cas(k2)-invtau_v_nudg_prof_cas(k1))
    1060 
    1061          w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
    1062          omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1))
    1063          du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
    1064          hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
    1065          vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
    1066          dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
    1067          hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
    1068          vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
    1069          dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
    1070          ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
    1071          vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
    1072          dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1))
    1073          hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1))
    1074          vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1))
    1075          dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
    1076          hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
    1077          vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
    1078          dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
    1079      
    1080          else !play>plev_prof_cas(1)
    1081 
    1082          k1=1
    1083          k2=2
    1084          print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2)
    1085          frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
    1086          frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
    1087          t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
    1088          theta_mod_cas(l)= frac1*th_prof_cas(k1) - frac2*th_prof_cas(k2)
    1089          if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
    1090          thv_mod_cas(l)= frac1*thv_prof_cas(k1) - frac2*thv_prof_cas(k2)
    1091          thl_mod_cas(l)= frac1*thl_prof_cas(k1) - frac2*thl_prof_cas(k2)
    1092          qv_mod_cas(l)= frac1*qv_prof_cas(k1) - frac2*qv_prof_cas(k2)
    1093          ql_mod_cas(l)= frac1*ql_prof_cas(k1) - frac2*ql_prof_cas(k2)
    1094          qi_mod_cas(l)= frac1*qi_prof_cas(k1) - frac2*qi_prof_cas(k2)
    1095          u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
    1096          v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
    1097          ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
    1098          vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
    1099          temp_nudg_mod_cas(l)= frac1*temp_nudg_prof_cas(k1) - frac2*temp_nudg_prof_cas(k2)
    1100          qv_nudg_mod_cas(l)= frac1*qv_nudg_prof_cas(k1) - frac2*qv_nudg_prof_cas(k2)
    1101          u_nudg_mod_cas(l)= frac1*u_nudg_prof_cas(k1) - frac2*u_nudg_prof_cas(k2)
    1102          v_nudg_mod_cas(l)= frac1*v_nudg_prof_cas(k1) - frac2*v_nudg_prof_cas(k2)
    1103 
    1104          invtau_temp_nudg_mod_cas(l)= frac1*invtau_temp_nudg_prof_cas(k1) - frac2*invtau_temp_nudg_prof_cas(k2)
    1105          invtau_qv_nudg_mod_cas(l)= frac1*invtau_qv_nudg_prof_cas(k1) - frac2*invtau_qv_nudg_prof_cas(k2)
    1106          invtau_u_nudg_mod_cas(l)= frac1*invtau_u_nudg_prof_cas(k1) - frac2*invtau_u_nudg_prof_cas(k2)
    1107          invtau_v_nudg_mod_cas(l)= frac1*invtau_v_nudg_prof_cas(k1) - frac2*invtau_v_nudg_prof_cas(k2)
    1108 
    1109          w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
    1110          omega_mod_cas(l)= frac1*omega_prof_cas(k1) - frac2*omega_prof_cas(k2)
    1111          du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
    1112          hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
    1113          vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
    1114          dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
    1115          hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
    1116          vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
    1117          dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
    1118          ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
    1119          vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
    1120          dth_mod_cas(l)= frac1*dth_prof_cas(k1) - frac2*dth_prof_cas(k2)
    1121          hth_mod_cas(l)= frac1*hth_prof_cas(k1) - frac2*hth_prof_cas(k2)
    1122          vth_mod_cas(l)= frac1*vth_prof_cas(k1) - frac2*vth_prof_cas(k2)
    1123          dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
    1124          hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
    1125          vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
    1126          dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
    1127 
    1128          endif ! play.le.plev_prof_cas(1)
    1129 
    1130         else ! above max altitude of forcing file
    1131  
    1132 !jyg
    1133          fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
    1134          fact = max(fact,0.)                                           !jyg
    1135          fact = exp(-fact)                                             !jyg
    1136          t_mod_cas(l)= t_prof_cas(nlev_cas)                            !jyg
    1137          theta_mod_cas(l)= th_prof_cas(nlev_cas)                       !jyg
    1138          thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
    1139          thl_mod_cas(l)= thl_prof_cas(nlev_cas)                        !jyg
    1140          qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact                     !jyg
    1141          ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact                     !jyg
    1142          qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact                     !jyg
    1143          u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                       !jyg
    1144          v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                       !jyg
    1145          ug_mod_cas(l)= ug_prof_cas(nlev_cas)                          !jyg
    1146          vg_mod_cas(l)= vg_prof_cas(nlev_cas)                          !jyg
    1147          temp_nudg_mod_cas(l)= temp_nudg_prof_cas(nlev_cas)            !jyg
    1148          qv_nudg_mod_cas(l)= qv_nudg_prof_cas(nlev_cas)                !jyg
    1149          u_nudg_mod_cas(l)= u_nudg_prof_cas(nlev_cas)                  !jyg
    1150          v_nudg_mod_cas(l)= v_nudg_prof_cas(nlev_cas)                  !jyg
    1151 
    1152          invtau_temp_nudg_mod_cas(l)= invtau_temp_nudg_prof_cas(nlev_cas)            !jyg
    1153          invtau_qv_nudg_mod_cas(l)= invtau_qv_nudg_prof_cas(nlev_cas)                !jyg
    1154          invtau_u_nudg_mod_cas(l)= invtau_u_nudg_prof_cas(nlev_cas)                  !jyg
    1155          invtau_v_nudg_mod_cas(l)= invtau_v_nudg_prof_cas(nlev_cas)                  !jyg
    1156 
    1157          thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
    1158          w_mod_cas(l)= 0.0                                             !jyg
    1159          omega_mod_cas(l)= 0.0                                         !jyg
    1160          du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
    1161          hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                     !jyg
    1162          vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                     !jyg
    1163          dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
    1164          hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                     !jyg
    1165          vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                     !jyg
    1166          dt_mod_cas(l)= dt_prof_cas(nlev_cas)
    1167          ht_mod_cas(l)= ht_prof_cas(nlev_cas)                          !jyg
    1168          vt_mod_cas(l)= vt_prof_cas(nlev_cas)                          !jyg
    1169          dth_mod_cas(l)= dth_prof_cas(nlev_cas)
    1170          hth_mod_cas(l)= hth_prof_cas(nlev_cas)                        !jyg
    1171          vth_mod_cas(l)= vth_prof_cas(nlev_cas)                        !jyg
    1172          dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
    1173          hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
    1174          vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
    1175          dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact               !jyg
    1176  
    1177         endif ! play
    1178  
    1179        enddo ! l
    1180 
    1181 ! for variables defined at layer interfaces (EV):
    1182 
    1183 
    1184        do l = 1, llm+1
    1185 
    1186         if (plev(l).ge.plev_prof_cas(nlev_cas)) then
    1187 
    1188          mxcalc=l
    1189          k1=0
    1190          k2=0
    1191 
    1192          if (plev(l).le.plev_prof_cas(1)) then
    1193 
    1194          do k = 1, nlev_cas-1
    1195           if (plev(l).le.plev_prof_cas(k).and. plev(l).gt.plev_prof_cas(k+1)) then
    1196             k1=k
    1197             k2=k+1
    1198           endif
    1199          enddo
    1200 
    1201          if (k1.eq.0 .or. k2.eq.0) then
    1202           write(*,*) 'PB! k1, k2 = ',k1,k2
    1203           write(*,*) 'l,plev(l) = ',l,plev(l)/100
    1204          do k = 1, nlev_cas-1
    1205           write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
    1206          enddo
    1207          endif
    1208 
    1209          frac = (plev_prof_cas(k2)-plev(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
    1210          tke_mod_cas(l)= tke_prof_cas(k2) - frac*(tke_prof_cas(k2)-tke_prof_cas(k1))
    1211          else !play>plev_prof_cas(1)
    1212          k1=1
    1213          k2=2
    1214          tke_mod_cas(l)= frac1*tke_prof_cas(k1) - frac2*tke_prof_cas(k2)
    1215 
    1216          endif ! plev.le.plev_prof_cas(1)
    1217 
    1218         else ! above max altitude of forcing file
    1219 
    1220          tke_mod_cas(l)=0.0
    1221 
    1222         endif ! plev
    1223 
    1224        enddo ! l
    1225 
    1226 
    1227 
    1228           return
    1229           end SUBROUTINE interp2_case_vertical_std
    1230 !*****************************************************************************
     939    !-------------------------------------------------------------------------
     940    ! Vertical interpolation of generic case forcing data onto mod_casel levels
     941    !-------------------------------------------------------------------------
     942
     943    integer nlevmax
     944    parameter (nlevmax=41)
     945    integer nlev_cas,mxcalc
     946    !       real play(llm), plev_prof(nlevmax)
     947    !       real t_prof(nlevmax),q_prof(nlevmax)
     948    !       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
     949    !       real ht_prof(nlevmax),vt_prof(nlevmax)
     950    !       real hq_prof(nlevmax),vq_prof(nlevmax)
     951
     952    real play(llm), plev(llm+1), plev_prof_cas(nlev_cas)
     953    real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
     954    real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
     955    real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
     956    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)
     957    real temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas)
     958    real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
     959    real invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas)
     960    real invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas)
     961
     962    real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
     963    real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
     964    real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
     965    real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
     966    real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
     967
     968    real t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm)
     969    real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
     970    real u_mod_cas(llm),v_mod_cas(llm)
     971    real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm),tke_mod_cas(llm+1)
     972    real temp_nudg_mod_cas(llm),qv_nudg_mod_cas(llm)
     973    real u_nudg_mod_cas(llm),v_nudg_mod_cas(llm)
     974    real invtau_temp_nudg_mod_cas(llm),invtau_qv_nudg_mod_cas(llm)
     975    real invtau_u_nudg_mod_cas(llm),invtau_v_nudg_mod_cas(llm)
     976    real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
     977    real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
     978    real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
     979    real dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm)
     980    real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
     981
     982    integer l,k,k1,k2
     983    real frac,frac1,frac2,fact
     984
     985
     986
     987    ! for variables defined at the middle of layers
     988
     989    do l = 1, llm
     990
     991       if (play(l).ge.plev_prof_cas(nlev_cas)) then
     992
     993          mxcalc=l
     994          !        print *,'debut interp2, mxcalc=',mxcalc
     995          k1=0
     996          k2=0
     997
     998          if (play(l).le.plev_prof_cas(1)) then
     999
     1000             do k = 1, nlev_cas-1
     1001                if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
     1002                   k1=k
     1003                   k2=k+1
     1004                endif
     1005             enddo
     1006
     1007             if (k1.eq.0 .or. k2.eq.0) then
     1008                write(*,*) 'PB! k1, k2 = ',k1,k2
     1009                write(*,*) 'l,play(l) = ',l,play(l)/100
     1010                do k = 1, nlev_cas-1
     1011                   write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
     1012                enddo
     1013             endif
     1014
     1015
     1016
     1017             frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
     1018
     1019             t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
     1020             theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
     1021             if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
     1022             thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1))
     1023             thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1))
     1024             qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1))
     1025             ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1))
     1026             qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1))
     1027             u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
     1028             v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
     1029             ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
     1030             vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
     1031             temp_nudg_mod_cas(l)= temp_nudg_prof_cas(k2) - frac*(temp_nudg_prof_cas(k2)-temp_nudg_prof_cas(k1))
     1032             qv_nudg_mod_cas(l)= qv_nudg_prof_cas(k2) - frac*(qv_nudg_prof_cas(k2)-qv_nudg_prof_cas(k1))
     1033             u_nudg_mod_cas(l)= u_nudg_prof_cas(k2) - frac*(u_nudg_prof_cas(k2)-u_nudg_prof_cas(k1))
     1034             v_nudg_mod_cas(l)= v_nudg_prof_cas(k2) - frac*(v_nudg_prof_cas(k2)-v_nudg_prof_cas(k1))
     1035
     1036             invtau_temp_nudg_mod_cas(l)= invtau_temp_nudg_prof_cas(k2) &
     1037                  - frac*(invtau_temp_nudg_prof_cas(k2)-invtau_temp_nudg_prof_cas(k1))
     1038             invtau_qv_nudg_mod_cas(l)= invtau_qv_nudg_prof_cas(k2) - frac*(invtau_qv_nudg_prof_cas(k2)-invtau_qv_nudg_prof_cas(k1))
     1039             invtau_u_nudg_mod_cas(l)= invtau_u_nudg_prof_cas(k2) - frac*(invtau_u_nudg_prof_cas(k2)-invtau_u_nudg_prof_cas(k1))
     1040             invtau_v_nudg_mod_cas(l)= invtau_v_nudg_prof_cas(k2) - frac*(invtau_v_nudg_prof_cas(k2)-invtau_v_nudg_prof_cas(k1))
     1041
     1042             w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
     1043             omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1))
     1044             du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
     1045             hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
     1046             vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
     1047             dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
     1048             hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
     1049             vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
     1050             dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
     1051             ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
     1052             vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
     1053             dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1))
     1054             hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1))
     1055             vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1))
     1056             dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
     1057             hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
     1058             vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
     1059             dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
     1060
     1061          else !play>plev_prof_cas(1)
     1062
     1063             k1=1
     1064             k2=2
     1065             print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2)
     1066             frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
     1067             frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
     1068             t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
     1069             theta_mod_cas(l)= frac1*th_prof_cas(k1) - frac2*th_prof_cas(k2)
     1070             if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
     1071             thv_mod_cas(l)= frac1*thv_prof_cas(k1) - frac2*thv_prof_cas(k2)
     1072             thl_mod_cas(l)= frac1*thl_prof_cas(k1) - frac2*thl_prof_cas(k2)
     1073             qv_mod_cas(l)= frac1*qv_prof_cas(k1) - frac2*qv_prof_cas(k2)
     1074             ql_mod_cas(l)= frac1*ql_prof_cas(k1) - frac2*ql_prof_cas(k2)
     1075             qi_mod_cas(l)= frac1*qi_prof_cas(k1) - frac2*qi_prof_cas(k2)
     1076             u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
     1077             v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
     1078             ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
     1079             vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
     1080             temp_nudg_mod_cas(l)= frac1*temp_nudg_prof_cas(k1) - frac2*temp_nudg_prof_cas(k2)
     1081             qv_nudg_mod_cas(l)= frac1*qv_nudg_prof_cas(k1) - frac2*qv_nudg_prof_cas(k2)
     1082             u_nudg_mod_cas(l)= frac1*u_nudg_prof_cas(k1) - frac2*u_nudg_prof_cas(k2)
     1083             v_nudg_mod_cas(l)= frac1*v_nudg_prof_cas(k1) - frac2*v_nudg_prof_cas(k2)
     1084
     1085             invtau_temp_nudg_mod_cas(l)= frac1*invtau_temp_nudg_prof_cas(k1) - frac2*invtau_temp_nudg_prof_cas(k2)
     1086             invtau_qv_nudg_mod_cas(l)= frac1*invtau_qv_nudg_prof_cas(k1) - frac2*invtau_qv_nudg_prof_cas(k2)
     1087             invtau_u_nudg_mod_cas(l)= frac1*invtau_u_nudg_prof_cas(k1) - frac2*invtau_u_nudg_prof_cas(k2)
     1088             invtau_v_nudg_mod_cas(l)= frac1*invtau_v_nudg_prof_cas(k1) - frac2*invtau_v_nudg_prof_cas(k2)
     1089
     1090             w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
     1091             omega_mod_cas(l)= frac1*omega_prof_cas(k1) - frac2*omega_prof_cas(k2)
     1092             du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
     1093             hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
     1094             vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
     1095             dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
     1096             hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
     1097             vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
     1098             dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
     1099             ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
     1100             vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
     1101             dth_mod_cas(l)= frac1*dth_prof_cas(k1) - frac2*dth_prof_cas(k2)
     1102             hth_mod_cas(l)= frac1*hth_prof_cas(k1) - frac2*hth_prof_cas(k2)
     1103             vth_mod_cas(l)= frac1*vth_prof_cas(k1) - frac2*vth_prof_cas(k2)
     1104             dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
     1105             hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
     1106             vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
     1107             dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
     1108
     1109          endif ! play.le.plev_prof_cas(1)
     1110
     1111       else ! above max altitude of forcing file
     1112
     1113          !jyg
     1114          fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
     1115          fact = max(fact,0.)                                           !jyg
     1116          fact = exp(-fact)                                             !jyg
     1117          t_mod_cas(l)= t_prof_cas(nlev_cas)                            !jyg
     1118          theta_mod_cas(l)= th_prof_cas(nlev_cas)                       !jyg
     1119          thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
     1120          thl_mod_cas(l)= thl_prof_cas(nlev_cas)                        !jyg
     1121          qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact                     !jyg
     1122          ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact                     !jyg
     1123          qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact                     !jyg
     1124          u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                       !jyg
     1125          v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                       !jyg
     1126          ug_mod_cas(l)= ug_prof_cas(nlev_cas)                          !jyg
     1127          vg_mod_cas(l)= vg_prof_cas(nlev_cas)                          !jyg
     1128          temp_nudg_mod_cas(l)= temp_nudg_prof_cas(nlev_cas)            !jyg
     1129          qv_nudg_mod_cas(l)= qv_nudg_prof_cas(nlev_cas)                !jyg
     1130          u_nudg_mod_cas(l)= u_nudg_prof_cas(nlev_cas)                  !jyg
     1131          v_nudg_mod_cas(l)= v_nudg_prof_cas(nlev_cas)                  !jyg
     1132
     1133          invtau_temp_nudg_mod_cas(l)= invtau_temp_nudg_prof_cas(nlev_cas)            !jyg
     1134          invtau_qv_nudg_mod_cas(l)= invtau_qv_nudg_prof_cas(nlev_cas)                !jyg
     1135          invtau_u_nudg_mod_cas(l)= invtau_u_nudg_prof_cas(nlev_cas)                  !jyg
     1136          invtau_v_nudg_mod_cas(l)= invtau_v_nudg_prof_cas(nlev_cas)                  !jyg
     1137
     1138          thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
     1139          w_mod_cas(l)= 0.0                                             !jyg
     1140          omega_mod_cas(l)= 0.0                                         !jyg
     1141          du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
     1142          hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                     !jyg
     1143          vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                     !jyg
     1144          dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
     1145          hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                     !jyg
     1146          vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                     !jyg
     1147          dt_mod_cas(l)= dt_prof_cas(nlev_cas)
     1148          ht_mod_cas(l)= ht_prof_cas(nlev_cas)                          !jyg
     1149          vt_mod_cas(l)= vt_prof_cas(nlev_cas)                          !jyg
     1150          dth_mod_cas(l)= dth_prof_cas(nlev_cas)
     1151          hth_mod_cas(l)= hth_prof_cas(nlev_cas)                        !jyg
     1152          vth_mod_cas(l)= vth_prof_cas(nlev_cas)                        !jyg
     1153          dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
     1154          hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
     1155          vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
     1156          dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact               !jyg
     1157
     1158       endif ! play
     1159
     1160    enddo ! l
     1161
     1162    ! for variables defined at layer interfaces (EV):
     1163
     1164
     1165    do l = 1, llm+1
     1166
     1167       if (plev(l).ge.plev_prof_cas(nlev_cas)) then
     1168
     1169          mxcalc=l
     1170          k1=0
     1171          k2=0
     1172
     1173          if (plev(l).le.plev_prof_cas(1)) then
     1174
     1175             do k = 1, nlev_cas-1
     1176                if (plev(l).le.plev_prof_cas(k).and. plev(l).gt.plev_prof_cas(k+1)) then
     1177                   k1=k
     1178                   k2=k+1
     1179                endif
     1180             enddo
     1181
     1182             if (k1.eq.0 .or. k2.eq.0) then
     1183                write(*,*) 'PB! k1, k2 = ',k1,k2
     1184                write(*,*) 'l,plev(l) = ',l,plev(l)/100
     1185                do k = 1, nlev_cas-1
     1186                   write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
     1187                enddo
     1188             endif
     1189
     1190             frac = (plev_prof_cas(k2)-plev(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
     1191             tke_mod_cas(l)= tke_prof_cas(k2) - frac*(tke_prof_cas(k2)-tke_prof_cas(k1))
     1192          else !play>plev_prof_cas(1)
     1193             k1=1
     1194             k2=2
     1195             frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
     1196             frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
     1197             tke_mod_cas(l)= frac1*tke_prof_cas(k1) - frac2*tke_prof_cas(k2)
     1198
     1199          endif ! plev.le.plev_prof_cas(1)
     1200
     1201       else ! above max altitude of forcing file
     1202
     1203          tke_mod_cas(l)=0.0
     1204
     1205       endif ! plev
     1206
     1207    enddo ! l
     1208
     1209
     1210
     1211    return
     1212  end SUBROUTINE interp2_case_vertical_std
     1213  !*****************************************************************************
    12311214
    12321215
Note: See TracChangeset for help on using the changeset viewer.