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

Merge revisions 3427:3600 of trunk into branch Ocean_skin

Location:
LMDZ6/branches/Ocean_skin
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Ocean_skin

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

    • Property svn:keywords set to Id
    r3316 r3605  
     1!
     2! $Id$
     3!
    14!#ifdef CPP_1D
    25!#include "../dyn3d/mod_const_mpi.F90"
    … …  
    69
    710
    8       PROGRAM lmdz1d
     11   PROGRAM lmdz1d
    912
    10    USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar
    11    USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
    12        clwcon, detr_therm, &
    13        qsol, fevap, z0m, z0h, agesno, &
    14        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
    15        falb_dir, falb_dif, &
    16        ftsol, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
    17        rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
    18        solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, &
    19        wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
    20        wake_deltaq, wake_deltat, wake_s, wake_dens, &
    21        zgam, zmax0, zmea, zpic, zsig, &
    22        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, ql_ancien, qs_ancien, &
    23        prlw_ancien, prsw_ancien, prw_ancien
    24  
    25    USE dimphy
    26    USE surface_data, only : type_ocean,ok_veget
    27    USE pbl_surface_mod, only : ftsoil, pbl_surface_init, &
    28                                  pbl_surface_final
    29    USE fonte_neige_mod, only : fonte_neige_init, fonte_neige_final
     13   USE ioipsl, only: getin
    3014
    31    USE infotrac ! new
    32    USE control_mod
    33    USE indice_sol_mod
    34    USE phyaqua_mod
    35 !  USE mod_1D_cases_read
    36    USE mod_1D_cases_read2
    37    USE mod_1D_amma_read
    38    USE print_control_mod, ONLY: lunout, prt_level
    39    USE iniphysiq_mod, ONLY: iniphysiq
    40    USE mod_const_mpi, ONLY: comm_lmdz
    41    USE physiq_mod, ONLY: physiq
    42    USE comvert_mod, ONLY: presnivs, ap, bp, dpres,nivsig, nivsigs, pa, &
    43                           preff, aps, bps, pseudoalt, scaleheight
    44    USE temps_mod, ONLY: annee_ref, calend, day_end, day_ini, day_ref, &
    45                         itau_dyn, itau_phy, start_time
     15   INTEGER forcing_type
    4616
    47       implicit none
    48 #include "dimensions.h"
    49 #include "YOMCST.h"
    50 !!#include "control.h"
    51 #include "clesphys.h"
    52 #include "dimsoil.h"
    53 !#include "indicesol.h"
     17   CALL getin('forcing_type',forcing_type)
    5418
    55 #include "compar1d.h"
    56 #include "flux_arp.h"
    57 #include "date_cas.h"
    58 #include "tsoilnudge.h"
    59 #include "fcg_gcssold.h"
    60 !!!#include "fbforcing.h"
    61 #include "compbl.h"
     19   IF (forcing_type==1000) THEN
     20           CALL scm
     21   ELSE
     22           CALL old_lmdz1d
     23   ENDIF
    6224
    63 !=====================================================================
    64 ! DECLARATIONS
    65 !=====================================================================
     25   END
    6626
    67 !---------------------------------------------------------------------
    68 !  Externals
    69 !---------------------------------------------------------------------
    70       external fq_sat
    71       real fq_sat
    72 
    73 !---------------------------------------------------------------------
    74 !  Arguments d' initialisations de la physique (USER DEFINE)
    75 !---------------------------------------------------------------------
    76 
    77       integer, parameter :: ngrid=1
    78       real :: zcufi    = 1.
    79       real :: zcvfi    = 1.
    80 
    81 !-      real :: nat_surf
    82 !-      logical :: ok_flux_surf
    83 !-      real :: fsens
    84 !-      real :: flat
    85 !-      real :: tsurf
    86 !-      real :: rugos
    87 !-      real :: qsol(1:2)
    88 !-      real :: qsurf
    89 !-      real :: psurf
    90 !-      real :: zsurf
    91 !-      real :: albedo
    92 !-
    93 !-      real :: time     = 0.
    94 !-      real :: time_ini
    95 !-      real :: xlat
    96 !-      real :: xlon
    97 !-      real :: wtsurf
    98 !-      real :: wqsurf
    99 !-      real :: restart_runoff
    100 !-      real :: xagesno
    101 !-      real :: qsolinp
    102 !-      real :: zpicinp
    103 !-
    104       real :: fnday
    105       real :: day, daytime
    106       real :: day1
    107       real :: heure
    108       integer :: jour
    109       integer :: mois
    110       integer :: an
    111  
    112 !---------------------------------------------------------------------
    113 !  Declarations related to forcing and initial profiles
    114 !---------------------------------------------------------------------
    115 
    116         integer :: kmax = llm
    117         integer llm700,nq1,nq2
    118         INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
    119         real timestep, frac
    120         real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max)
    121         real  uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max)
    122         real  ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max)
    123         real  dqtdxls(nlev_max),dqtdyls(nlev_max)
    124         real  dqtdtls(nlev_max),thlpcar(nlev_max)
    125         real  qprof(nlev_max,nqmx)
    126 
    127 !        integer :: forcing_type
    128         logical :: forcing_les     = .false.
    129         logical :: forcing_armcu   = .false.
    130         logical :: forcing_rico    = .false.
    131         logical :: forcing_radconv = .false.
    132         logical :: forcing_toga    = .false.
    133         logical :: forcing_twpice  = .false.
    134         logical :: forcing_amma    = .false.
    135         logical :: forcing_dice    = .false.
    136         logical :: forcing_gabls4  = .false.
    137 
    138         logical :: forcing_GCM2SCM = .false.
    139         logical :: forcing_GCSSold = .false.
    140         logical :: forcing_sandu   = .false.
    141         logical :: forcing_astex   = .false.
    142         logical :: forcing_fire    = .false.
    143         logical :: forcing_case    = .false.
    144         logical :: forcing_case2   = .false.
    145         integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
    146 !                                                            (cf read_tsurf1d.F)
    147 
    148 !vertical advection computation
    149 !       real d_t_z(llm), d_q_z(llm)
    150 !       real d_t_dyn_z(llm), dq_dyn_z(llm)
    151 !       real zz(llm)
    152 !       real zfact
    153 
    154 !flag forcings
    155         logical :: nudge_wind=.true.
    156         logical :: nudge_thermo=.false.
    157         logical :: cptadvw=.true.
    158 !=====================================================================
    159 ! DECLARATIONS FOR EACH CASE
    160 !=====================================================================
    161 !
    162 #include "1D_decl_cases.h"
    163 !
    164 !---------------------------------------------------------------------
    165 !  Declarations related to nudging
    166 !---------------------------------------------------------------------
    167      integer :: nudge_max
    168      parameter (nudge_max=9)
    169      integer :: inudge_RHT=1
    170      integer :: inudge_UV=2
    171      logical :: nudge(nudge_max)
    172      real :: t_targ(llm)
    173      real :: rh_targ(llm)
    174      real :: u_targ(llm)
    175      real :: v_targ(llm)
    176 !
    177 !---------------------------------------------------------------------
    178 !  Declarations related to vertical discretization:
    179 !---------------------------------------------------------------------
    180       real :: pzero=1.e5
    181       real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
    182       real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1)
    183 
    184 !---------------------------------------------------------------------
    185 !  Declarations related to variables
    186 !---------------------------------------------------------------------
    187 
    188       real :: phi(llm)
    189       real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
    190       REAL rot(1, llm) ! relative vorticity, in s-1
    191       real :: rlat_rad(1),rlon_rad(1)
    192       real :: omega(llm+1),omega2(llm),rho(llm+1)
    193       real :: ug(llm),vg(llm),fcoriolis
    194       real :: sfdt, cfdt
    195       real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
    196       real :: dt_dyn(llm)
    197       real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
    198       real :: d_u_nudge(llm),d_v_nudge(llm)
    199       real :: du_adv(llm),dv_adv(llm)
    200       real :: du_age(llm),dv_age(llm)
    201       real :: alpha
    202       real :: ttt
    203 
    204       REAL, ALLOCATABLE, DIMENSION(:,:):: q
    205       REAL, ALLOCATABLE, DIMENSION(:,:):: dq
    206       REAL, ALLOCATABLE, DIMENSION(:,:):: dq_dyn
    207       REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
    208       REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
    209 !      REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
    210 
    211 !---------------------------------------------------------------------
    212 !  Initialization of surface variables
    213 !---------------------------------------------------------------------
    214       real :: run_off_lic_0(1)
    215       real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
    216       real :: tsoil(1,nsoilmx,nbsrf)
    217 !     real :: agesno(1,nbsrf)
    218 
    219 !---------------------------------------------------------------------
    220 !  Call to phyredem
    221 !---------------------------------------------------------------------
    222       logical :: ok_writedem =.true.
    223       real :: sollw_in = 0.
    224       real :: solsw_in = 0.
    225      
    226 !---------------------------------------------------------------------
    227 !  Call to physiq
    228 !---------------------------------------------------------------------
    229       logical :: firstcall=.true.
    230       logical :: lastcall=.false.
    231       real :: phis(1)    = 0.0
    232       real :: dpsrf(1)
    233 
    234 !---------------------------------------------------------------------
    235 !  Initializations of boundary conditions
    236 !---------------------------------------------------------------------
    237       integer, parameter :: yd = 360
    238       real :: phy_nat (yd) = 0.0 ! 0=ocean libre,1=land,2=glacier,3=banquise
    239       real :: phy_alb (yd)  ! Albedo land only (old value condsurf_jyg=0.3)
    240       real :: phy_sst (yd)  ! SST (will not be used; cf read_tsurf1d.F)
    241       real :: phy_bil (yd) = 1.0 ! Ne sert que pour les slab_ocean
    242       real :: phy_rug (yd) ! Longueur rugosite utilisee sur land only
    243       real :: phy_ice (yd) = 0.0 ! Fraction de glace
    244       real :: phy_fter(yd) = 0.0 ! Fraction de terre
    245       real :: phy_foce(yd) = 0.0 ! Fraction de ocean
    246       real :: phy_fsic(yd) = 0.0 ! Fraction de glace
    247       real :: phy_flic(yd) = 0.0 ! Fraction de glace
    248 
    249 !---------------------------------------------------------------------
    250 !  Fichiers et d'autres variables
    251 !---------------------------------------------------------------------
    252       integer :: k,l,i,it=1,mxcalc
    253       integer :: nsrf
    254       integer jcode
    255       INTEGER read_climoz
    256 !
    257       integer :: it_end ! iteration number of the last call
    258 !Al1
    259       integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
    260       data ecrit_slab_oc/-1/
    261 !
    262 !     if flag_inhib_forcing = 0, tendencies of forcing are added
    263 !                           <> 0, tendencies of forcing are not added
    264       INTEGER :: flag_inhib_forcing = 0
    265 
    266 !=====================================================================
    267 ! INITIALIZATIONS
    268 !=====================================================================
    269       du_phys(:)=0.
    270       dv_phys(:)=0.
    271       dt_phys(:)=0.
    272       dt_dyn(:)=0.
    273       dt_cooling(:)=0.
    274       d_t_adv(:)=0.
    275       d_t_nudge(:)=0.
    276       d_u_nudge(:)=0.
    277       d_v_nudge(:)=0.
    278       du_adv(:)=0.
    279       dv_adv(:)=0.
    280       du_age(:)=0.
    281       dv_age(:)=0.
    282      
    283 ! Initialization of Common turb_forcing
    284        dtime_frcg = 0.
    285        Turb_fcg_gcssold=.false.
    286        hthturb_gcssold = 0.
    287        hqturb_gcssold = 0.
    288 
    289 
    290 
    291 
    292 !---------------------------------------------------------------------
    293 ! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
    294 !---------------------------------------------------------------------
    295 !Al1
    296         call conf_unicol
    297 !Al1 moves this gcssold var from common fcg_gcssold to
    298         Turb_fcg_gcssold = xTurb_fcg_gcssold
    299 ! --------------------------------------------------------------------
    300         close(1)
    301 !Al1
    302         write(*,*) 'lmdz1d.def lu => unicol.def'
    303 
    304 ! forcing_type defines the way the SCM is forced:
    305 !forcing_type = 0 ==> forcing_les = .true.
    306 !             initial profiles from file prof.inp.001
    307 !             no forcing by LS convergence ;
    308 !             surface temperature imposed ;
    309 !             radiative cooling may be imposed (iflag_radia=0 in physiq.def)
    310 !forcing_type = 1 ==> forcing_radconv = .true.
    311 !             idem forcing_type = 0, but the imposed radiative cooling
    312 !             is set to 0 (hence, if iflag_radia=0 in physiq.def,
    313 !             then there is no radiative cooling at all)
    314 !forcing_type = 2 ==> forcing_toga = .true.
    315 !             initial profiles from TOGA-COARE IFA files
    316 !             LS convergence and SST imposed from TOGA-COARE IFA files
    317 !forcing_type = 3 ==> forcing_GCM2SCM = .true.
    318 !             initial profiles from the GCM output
    319 !             LS convergence imposed from the GCM output
    320 !forcing_type = 4 ==> forcing_twpice = .true.
    321 !             initial profiles from TWP-ICE cdf file
    322 !             LS convergence, omega and SST imposed from TWP-ICE files
    323 !forcing_type = 5 ==> forcing_rico = .true.
    324 !             initial profiles from RICO files
    325 !             LS convergence imposed from RICO files
    326 !forcing_type = 6 ==> forcing_amma = .true.
    327 !             initial profiles from AMMA nc file
    328 !             LS convergence, omega and surface fluxes imposed from AMMA file 
    329 !forcing_type = 7 ==> forcing_dice = .true.
    330 !             initial profiles and large scale forcings in dice_driver.nc
    331 !             Different stages: soil model alone, atm. model alone
    332 !             then both models coupled
    333 !forcing_type = 8 ==> forcing_gabls4 = .true.
    334 !             initial profiles and large scale forcings in gabls4_driver.nc
    335 !forcing_type >= 100 ==> forcing_case = .true.
    336 !             initial profiles and large scale forcings in cas.nc
    337 !             LS convergence, omega and SST imposed from CINDY-DYNAMO files
    338 !             101=cindynamo
    339 !             102=bomex
    340 !forcing_type >= 100 ==> forcing_case2 = .true.
    341 !             temporary flag while all the 1D cases are not whith the same cas.nc forcing file
    342 !             103=arm_cu2 ie arm_cu with new forcing format
    343 !             104=rico2 ie rico with new forcing format
    344 !forcing_type = 40 ==> forcing_GCSSold = .true.
    345 !             initial profile from GCSS file
    346 !             LS convergence imposed from GCSS file
    347 !forcing_type = 50 ==> forcing_fire = .true.
    348 !             forcing from fire.nc
    349 !forcing_type = 59 ==> forcing_sandu = .true.
    350 !             initial profiles from sanduref file: see prof.inp.001
    351 !             SST varying with time and divergence constante: see ifa_sanduref.txt file
    352 !             Radiation has to be computed interactively
    353 !forcing_type = 60 ==> forcing_astex = .true.
    354 !             initial profiles from file: see prof.inp.001
    355 !             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
    356 !             Radiation has to be computed interactively
    357 !forcing_type = 61 ==> forcing_armcu = .true.
    358 !             initial profiles from file: see prof.inp.001
    359 !             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
    360 !             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
    361 !             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
    362 !             Radiation to be switched off
    363 !
    364       if (forcing_type <=0) THEN
    365        forcing_les = .true.
    366       elseif (forcing_type .eq.1) THEN
    367        forcing_radconv = .true.
    368       elseif (forcing_type .eq.2) THEN
    369        forcing_toga    = .true.
    370       elseif (forcing_type .eq.3) THEN
    371        forcing_GCM2SCM = .true.
    372       elseif (forcing_type .eq.4) THEN
    373        forcing_twpice = .true.
    374       elseif (forcing_type .eq.5) THEN
    375        forcing_rico = .true.
    376       elseif (forcing_type .eq.6) THEN
    377        forcing_amma = .true.
    378       elseif (forcing_type .eq.7) THEN
    379        forcing_dice = .true.
    380       elseif (forcing_type .eq.8) THEN
    381        forcing_gabls4 = .true.
    382       elseif (forcing_type .eq.101) THEN ! Cindynamo starts 1-10-2011 0h
    383        forcing_case = .true.
    384        year_ini_cas=2011
    385        mth_ini_cas=10
    386        day_deb=1
    387        heure_ini_cas=0.
    388        pdt_cas=3*3600.         ! forcing frequency
    389       elseif (forcing_type .eq.102) THEN ! Bomex starts 24-6-1969 0h
    390        forcing_case = .true.
    391        year_ini_cas=1969
    392        mth_ini_cas=6
    393        day_deb=24
    394        heure_ini_cas=0.
    395        pdt_cas=1800.         ! forcing frequency
    396       elseif (forcing_type .eq.103) THEN ! Arm_cu starts 21-6-1997 11h30
    397        forcing_case2 = .true.
    398        year_ini_cas=1997
    399        mth_ini_cas=6
    400        day_deb=21
    401        heure_ini_cas=11.5
    402        pdt_cas=1800.         ! forcing frequency
    403       elseif (forcing_type .eq.104) THEN ! rico starts 16-12-2004 0h
    404        forcing_case2 = .true.
    405        year_ini_cas=2004
    406        mth_ini_cas=12
    407        day_deb=16
    408        heure_ini_cas=0.
    409        pdt_cas=1800.         ! forcing frequency
    410       elseif (forcing_type .eq.105) THEN ! bomex starts 16-12-2004 0h
    411        forcing_case2 = .true.
    412        year_ini_cas=1969
    413        mth_ini_cas=6
    414        day_deb=24
    415        heure_ini_cas=0.
    416        pdt_cas=1800.         ! forcing frequency
    417       elseif (forcing_type .eq.106) THEN ! ayotte_24SC starts 6-11-1992 0h
    418        forcing_case2 = .true.
    419        year_ini_cas=1992
    420        mth_ini_cas=11
    421        day_deb=6
    422        heure_ini_cas=10.
    423        pdt_cas=86400.        ! forcing frequency
    424       elseif (forcing_type .eq.40) THEN
    425        forcing_GCSSold = .true.
    426       elseif (forcing_type .eq.50) THEN
    427        forcing_fire = .true.
    428       elseif (forcing_type .eq.59) THEN
    429        forcing_sandu   = .true.
    430       elseif (forcing_type .eq.60) THEN
    431        forcing_astex   = .true.
    432       elseif (forcing_type .eq.61) THEN
    433        forcing_armcu = .true.
    434        IF(llm.NE.19.AND.llm.NE.40) stop 'Erreur nombre de niveaux !!'
    435       else
    436        write (*,*) 'ERROR : unknown forcing_type ', forcing_type
    437        stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61'
    438       ENDIF
    439       print*,"forcing type=",forcing_type
    440 
    441 ! if type_ts_forcing=0, the surface temp of 1D simulation is constant in time
    442 ! (specified by tsurf in lmdz1d.def); if type_ts_forcing=1, the surface temperature
    443 ! varies in time according to a forcing (e.g. forcing_toga) and is passed to read_tsurf1d.F
    444 ! through the common sst_forcing.
    445 
    446         type_ts_forcing = 0
    447         if (forcing_toga.or.forcing_sandu.or.forcing_astex .or. forcing_dice)                 &
    448      &    type_ts_forcing = 1
    449 !
    450 ! Initialization of the logical switch for nudging
    451      jcode = iflag_nudge
    452      do i = 1,nudge_max
    453        nudge(i) = mod(jcode,10) .ge. 1
    454        jcode = jcode/10
    455      enddo
    456 !---------------------------------------------------------------------
    457 !  Definition of the run
    458 !---------------------------------------------------------------------
    459 
    460       call conf_gcm( 99, .TRUE. )
    461 !-----------------------------------------------------------------------
    462 !   Choix du calendrier
    463 !   -------------------
    464 
    465 !      calend = 'earth_365d'
    466       if (calend == 'earth_360d') then
    467         call ioconf_calendar('360d')
    468         write(*,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
    469       else if (calend == 'earth_365d') then
    470         call ioconf_calendar('noleap')
    471         write(*,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
    472       else if (calend == 'earth_366d') then
    473         call ioconf_calendar('all_leap')
    474         write(*,*)'CALENDRIER CHOISI: Terrestre bissextile'
    475       else if (calend == 'gregorian') then
    476         call ioconf_calendar('gregorian') ! not to be used by normal users
    477         write(*,*)'CALENDRIER CHOISI: Gregorien'
    478       else
    479         write (*,*) 'ERROR : unknown calendar ', calend
    480         stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
    481       endif
    482 !-----------------------------------------------------------------------
    483 !
    484 !c Date :
    485 !      La date est supposee donnee sous la forme [annee, numero du jour dans
    486 !      l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
    487 !      On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
    488 !      Le numero du jour est dans "day". L heure est traitee separement.
    489 !      La date complete est dans "daytime" (l'unite est le jour).
    490       if (nday>0) then
    491          fnday=nday
    492       else
    493          fnday=-nday/float(day_step)
    494       endif
    495       print *,'fnday=',fnday
    496 !     start_time doit etre en FRACTION DE JOUR
    497       start_time=time_ini/24.
    498 
    499 ! Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
    500       IF(forcing_type .EQ. 61) fnday=53100./86400.
    501       IF(forcing_type .EQ. 103) fnday=53100./86400.
    502 ! Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
    503       IF(forcing_type .EQ. 6) fnday=64800./86400.
    504 !     IF(forcing_type .EQ. 6) fnday=50400./86400.
    505  IF(forcing_type .EQ. 8 ) fnday=129600./86400.
    506       annee_ref = anneeref
    507       mois = 1
    508       day_ref = dayref
    509       heure = 0.
    510       itau_dyn = 0
    511       itau_phy = 0
    512       call ymds2ju(annee_ref,mois,day_ref,heure,day)
    513       day_ini = int(day)
    514       day_end = day_ini + int(fnday)
    515 
    516       IF (forcing_type .eq.2) THEN
    517 ! Convert the initial date of Toga-Coare to Julian day
    518       call ymds2ju                                                          &
    519      & (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
    520 
    521       ELSEIF (forcing_type .eq.4) THEN
    522 ! Convert the initial date of TWPICE to Julian day
    523       call ymds2ju                                                          &
    524      & (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi              &
    525      & ,day_ju_ini_twpi)
    526       ELSEIF (forcing_type .eq.6) THEN
    527 ! Convert the initial date of AMMA to Julian day
    528       call ymds2ju                                                          &
    529      & (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma              &
    530      & ,day_ju_ini_amma)
    531       ELSEIF (forcing_type .eq.7) THEN
    532 ! Convert the initial date of DICE to Julian day
    533       call ymds2ju                                                         &
    534      & (year_ini_dice,mth_ini_dice,day_ini_dice,heure_ini_dice             &
    535      & ,day_ju_ini_dice)
    536  ELSEIF (forcing_type .eq.8 ) THEN
    537 ! Convert the initial date of GABLS4 to Julian day
    538       call ymds2ju                                                         &
    539      & (year_ini_gabls4,mth_ini_gabls4,day_ini_gabls4,heure_ini_gabls4     &
    540      & ,day_ju_ini_gabls4)
    541       ELSEIF (forcing_type .gt.100) THEN
    542 ! Convert the initial date to Julian day
    543       day_ini_cas=day_deb
    544       print*,'time case',year_ini_cas,mth_ini_cas,day_ini_cas
    545       call ymds2ju                                                         &
    546      & (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas*3600            &
    547      & ,day_ju_ini_cas)
    548       print*,'time case 2',day_ini_cas,day_ju_ini_cas
    549       ELSEIF (forcing_type .eq.59) THEN
    550 ! Convert the initial date of Sandu case to Julian day
    551       call ymds2ju                                                          &
    552      &   (year_ini_sandu,mth_ini_sandu,day_ini_sandu,                       &
    553      &    time_ini*3600.,day_ju_ini_sandu)
    554 
    555       ELSEIF (forcing_type .eq.60) THEN
    556 ! Convert the initial date of Astex case to Julian day
    557       call ymds2ju                                                          &
    558      &   (year_ini_astex,mth_ini_astex,day_ini_astex,                        &
    559      &    time_ini*3600.,day_ju_ini_astex)
    560 
    561       ELSEIF (forcing_type .eq.61) THEN
    562 ! Convert the initial date of Arm_cu case to Julian day
    563       call ymds2ju                                                          &
    564      & (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu          &
    565      & ,day_ju_ini_armcu)
    566       ENDIF
    567 
    568       IF (forcing_type .gt.100) THEN
    569       daytime = day + heure_ini_cas/24. ! 1st day and initial time of the simulation
    570       ELSE
    571       daytime = day + time_ini/24. ! 1st day and initial time of the simulation
    572       ENDIF
    573 ! Print out the actual date of the beginning of the simulation :
    574       call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
    575       print *,' Time of beginning : ',                                      &
    576      &        year_print, month_print, day_print, sec_print
    577 
    578 !---------------------------------------------------------------------
    579 ! Initialization of dimensions, geometry and initial state
    580 !---------------------------------------------------------------------
    581 !      call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
    582 !     but we still need to initialize dimphy module (klon,klev,etc.)  here.
    583       call init_dimphy(1,llm)
    584       call suphel
    585       call infotrac_init
    586 
    587       if (nqtot>nqmx) STOP 'Augmenter nqmx dans lmdz1d.F'
    588       allocate(q(llm,nqtot)) ; q(:,:)=0.
    589       allocate(dq(llm,nqtot))
    590       allocate(dq_dyn(llm,nqtot))
    591       allocate(d_q_adv(llm,nqtot))
    592       allocate(d_q_nudge(llm,nqtot))
    593 !      allocate(d_th_adv(llm))
    594 
    595       q(:,:) = 0.
    596       dq(:,:) = 0.
    597       dq_dyn(:,:) = 0.
    598       d_q_adv(:,:) = 0.
    599       d_q_nudge(:,:) = 0.
    600 
    601 !
    602 !   No ozone climatology need be read in this pre-initialization
    603 !          (phys_state_var_init is called again in physiq)
    604       read_climoz = 0
    605 !
    606       call phys_state_var_init(read_climoz)
    607 
    608       if (ngrid.ne.klon) then
    609          print*,'stop in inifis'
    610          print*,'Probleme de dimensions :'
    611          print*,'ngrid = ',ngrid
    612          print*,'klon  = ',klon
    613          stop
    614       endif
    615 !!!=====================================================================
    616 !!! Feedback forcing values for Gateaux differentiation (al1)
    617 !!!=====================================================================
    618 !!! Surface Planck forcing bracketing call radiation
    619 !!      surf_Planck = 0.
    620 !!      surf_Conv   = 0.
    621 !!      write(*,*) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
    622 !!! a mettre dans le lmdz1d.def ou autre
    623 !!
    624 !!
    625       qsol = qsolinp
    626       qsurf = fq_sat(tsurf,psurf/100.)
    627       day1= day_ini
    628       time=daytime-day
    629       ts_toga(1)=tsurf ! needed by read_tsurf1d.F
    630       rho(1)=psurf/(rd*tsurf*(1.+(rv/rd-1.)*qsurf))
    631 
    632 !
    633 !! mpl et jyg le 22/08/2012 :
    634 !!  pour que les cas a flux de surface imposes marchent
    635       IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
    636        fsens=-wtsurf*rcpd*rho(1)
    637        flat=-wqsurf*rlvtt*rho(1)
    638        print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
    639       ENDIF
    640       print*,'Flux sol ',fsens,flat
    641 !!      ok_flux_surf=.false.
    642 !!      fsens=-wtsurf*rcpd*rho(1)
    643 !!      flat=-wqsurf*rlvtt*rho(1)
    644 !!!!
    645 
    646 ! Vertical discretization and pressure levels at half and mid levels:
    647 
    648       pa   = 5e4
    649 !!      preff= 1.01325e5
    650       preff = psurf
    651       IF (ok_old_disvert) THEN
    652         call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    653         print *,'On utilise disvert0'
    654         aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
    655         bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
    656         scaleheight=8.
    657         pseudoalt(1:llm)=-scaleheight*log(presnivs(1:llm)/preff)
    658       ELSE
    659         call disvert()
    660         print *,'On utilise disvert'
    661 !       Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
    662 !       Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
    663       ENDIF
    664 
    665       sig_s=presnivs/preff
    666       plev =ap+bp*psurf
    667       play = 0.5*(plev(1:llm)+plev(2:llm+1))
    668       zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles
    669 
    670       IF (forcing_type .eq. 59) THEN
    671 ! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
    672       write(*,*) '***********************'
    673       do l = 1, llm
    674        write(*,*) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
    675        if (trouve_700 .and. play(l).le.70000) then
    676          llm700=l
    677          print *,'llm700,play=',llm700,play(l)/100.
    678          trouve_700= .false.
    679        endif
    680       enddo
    681       write(*,*) '***********************'
    682       ENDIF
    683 
    684 !
    685 !=====================================================================
    686 ! EVENTUALLY, READ FORCING DATA :
    687 !=====================================================================
    688 
    689 #include "1D_read_forc_cases.h"
    690 
    691       if (forcing_GCM2SCM) then
    692         write (*,*) 'forcing_GCM2SCM not yet implemented'
    693         stop 'in initialization'
    694       endif ! forcing_GCM2SCM
    695 
    696       print*,'mxcalc=',mxcalc
    697 !     print*,'zlay=',zlay(mxcalc)
    698       print*,'play=',play(mxcalc)
    699 
    700 !Al1 pour SST forced, appell?? depuis ocean_forced_noice
    701       ts_cur = tsurf ! SST used in read_tsurf1d
    702 !=====================================================================
    703 ! Initialisation de la physique :
    704 !=====================================================================
    705 
    706 !  Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
    707 !
    708 ! day_step, iphysiq lus dans gcm.def ci-dessus
    709 ! timestep: calcule ci-dessous from rday et day_step
    710 ! ngrid=1
    711 ! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
    712 ! rday: defini dans suphel.F (86400.)
    713 ! day_ini: lu dans run.def (dayref)
    714 ! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
    715 ! airefi,zcufi,zcvfi initialises au debut de ce programme
    716 ! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
    717       day_step = float(nsplit_phys)*day_step/float(iphysiq)
    718       write (*,*) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
    719       timestep =rday/day_step
    720       dtime_frcg = timestep
    721 !
    722       zcufi=airefi
    723       zcvfi=airefi
    724 !
    725       rlat_rad(1)=xlat*rpi/180.
    726       rlon_rad(1)=xlon*rpi/180.
    727 
    728      ! Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid,
    729      ! e.g. for cell boundaries, which are meaningless in 1D; so pad these
    730      ! with '0.' when necessary
    731       call iniphysiq(iim,jjm,llm, &
    732            1,comm_lmdz, &
    733            rday,day_ini,timestep,  &
    734            (/rlat_rad(1),0./),(/0./), &
    735            (/0.,0./),(/rlon_rad(1),0./),  &
    736            (/ (/airefi,0./),(/0.,0./) /), &
    737            (/zcufi,0.,0.,0./), &
    738            (/zcvfi,0./), &
    739            ra,rg,rd,rcpd,1)
    740       print*,'apres iniphysiq'
    741 
    742 ! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
    743       co2_ppm= 330.0
    744       solaire=1370.0
    745 
    746 ! Ecriture du startphy avant le premier appel a la physique.
    747 ! On le met juste avant pour avoir acces a tous les champs
    748 
    749       if (ok_writedem) then
    750 
    751 !--------------------------------------------------------------------------
    752 ! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
    753 ! need : qsol fder snow qsurf evap rugos agesno ftsoil
    754 !--------------------------------------------------------------------------
    755 
    756         type_ocean = "force"
    757         run_off_lic_0(1) = restart_runoff
    758         call fonte_neige_init(run_off_lic_0)
    759 
    760         fder=0.
    761         snsrf(1,:)=snowmass ! masse de neige des sous surface
    762         qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
    763         fevap=0.
    764         z0m(1,:)=rugos     ! couverture de neige des sous surface
    765         z0h(1,:)=rugosh    ! couverture de neige des sous surface
    766         agesno  = xagesno
    767         tsoil(:,:,:)=tsurf
    768 !------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
    769 !       tsoil(1,1,1)=299.18
    770 !       tsoil(1,2,1)=300.08
    771 !       tsoil(1,3,1)=301.88
    772 !       tsoil(1,4,1)=305.48
    773 !       tsoil(1,5,1)=308.00
    774 !       tsoil(1,6,1)=308.00
    775 !       tsoil(1,7,1)=308.00
    776 !       tsoil(1,8,1)=308.00
    777 !       tsoil(1,9,1)=308.00
    778 !       tsoil(1,10,1)=308.00
    779 !       tsoil(1,11,1)=308.00
    780 !-----------------------------------------------------------------------
    781         call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
    782 
    783 !------------------ prepare limit conditions for limit.nc -----------------
    784 !--   Ocean force
    785 
    786         print*,'avant phyredem'
    787         pctsrf(1,:)=0.
    788           if (nat_surf.eq.0.) then
    789           pctsrf(1,is_oce)=1.
    790           pctsrf(1,is_ter)=0.
    791           pctsrf(1,is_lic)=0.
    792           pctsrf(1,is_sic)=0.
    793         else if (nat_surf .eq. 1) then
    794           pctsrf(1,is_oce)=0.
    795           pctsrf(1,is_ter)=1.
    796           pctsrf(1,is_lic)=0.
    797           pctsrf(1,is_sic)=0.
    798         else if (nat_surf .eq. 2) then
    799           pctsrf(1,is_oce)=0.
    800           pctsrf(1,is_ter)=0.
    801           pctsrf(1,is_lic)=1.
    802           pctsrf(1,is_sic)=0.
    803         else if (nat_surf .eq. 3) then
    804           pctsrf(1,is_oce)=0.
    805           pctsrf(1,is_ter)=0.
    806           pctsrf(1,is_lic)=0.
    807           pctsrf(1,is_sic)=1.
    808 
    809      end if
    810 
    811 
    812         print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
    813      &        ,pctsrf(1,is_oce),pctsrf(1,is_ter)
    814 
    815         zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
    816         zpic = zpicinp
    817         ftsol=tsurf
    818         nsw=6 ! on met le nb de bandes SW=6, pour initialiser
    819               ! 6 albedo, mais on peut quand meme tourner avec
    820               ! moins. Seules les 2 ou 4 premiers seront lus
    821         falb_dir=albedo
    822         falb_dif=albedo
    823         rugoro=rugos
    824         t_ancien(1,:)=temp(:)
    825         q_ancien(1,:)=q(:,1)
    826         ql_ancien = 0.
    827         qs_ancien = 0.
    828         prlw_ancien = 0.
    829         prsw_ancien = 0.
    830         prw_ancien = 0.
    831 !jyg<
    832 !!        pbl_tke(:,:,:)=1.e-8
    833         pbl_tke(:,:,:)=0.
    834         pbl_tke(:,2,:)=1.e-2
    835         PRINT *, ' pbl_tke dans lmdz1d '
    836         if (prt_level .ge. 5) then
    837          DO nsrf = 1,4
    838            PRINT *,'pbl_tke(1,:,',nsrf,') ',pbl_tke(1,:,nsrf)
    839          ENDDO
    840         end if
    841 
    842 !>jyg
    843 
    844         rain_fall=0.
    845         snow_fall=0.
    846         solsw=0.
    847         sollw=0.
    848         sollwdown=rsigma*tsurf**4
    849         radsol=0.
    850         rnebcon=0.
    851         ratqs=0.
    852         clwcon=0.
    853         zmax0 = 0.
    854         zmea=0.
    855         zstd=0.
    856         zsig=0.
    857         zgam=0.
    858         zval=0.
    859         zthe=0.
    860         sig1=0.
    861         w01=0.
    862         wake_cstar = 0.
    863         wake_deltaq = 0.
    864         wake_deltat = 0.
    865         wake_delta_pbl_TKE(:,:,:) = 0.
    866         delta_tsurf = 0.
    867         wake_fip = 0.
    868         wake_pe = 0.
    869         wake_s = 0.
    870         wake_dens = 0.
    871         ale_bl = 0.
    872         ale_bl_trig = 0.
    873         alp_bl = 0.
    874         IF (ALLOCATED(du_gwd_rando)) du_gwd_rando = 0.
    875         IF (ALLOCATED(du_gwd_front)) du_gwd_front = 0.
    876         entr_therm = 0.
    877         detr_therm = 0.
    878         f0 = 0.
    879         fm_therm = 0.
    880         u_ancien(1,:)=u(:)
    881         v_ancien(1,:)=v(:)
    882  
    883 !------------------------------------------------------------------------
    884 ! Make file containing restart for the physics (startphy.nc)
    885 !
    886 ! NB: List of the variables to be written by phyredem (via put_field):
    887 ! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
    888 ! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
    889 ! qsol,falb_dir(:,nsrf),falb_dif(:,nsrf),evap(:,nsrf),snow(:,nsrf)
    890 ! radsol,solsw,sollw, sollwdown,fder,rain_fall,snow_fall,frugs(:,nsrf)
    891 ! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
    892 ! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
    893 ! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
    894 ! wake_deltat,wake_deltaq,wake_s,wake_dens,wake_cstar,
    895 ! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
    896 !
    897 ! NB2: The content of the startphy.nc file depends on some flags defined in
    898 ! the ".def" files. However, since conf_phys is not called in lmdz1d.F90, these flags have
    899 ! to be set at some arbitratry convenient values.
    900 !------------------------------------------------------------------------
    901 !Al1 =============== restart option ==========================
    902         if (.not.restart) then
    903           iflag_pbl = 5
    904           call phyredem ("startphy.nc")
    905         else
    906 ! (desallocations)
    907         print*,'callin surf final'
    908           call pbl_surface_final( fder, snsrf, qsurfsrf, tsoil)
    909         print*,'after surf final'
    910           CALL fonte_neige_final(run_off_lic_0)
    911         endif
    912 
    913         ok_writedem=.false.
    914         print*,'apres phyredem'
    915 
    916       endif ! ok_writedem
    917      
    918 !------------------------------------------------------------------------
    919 ! Make file containing boundary conditions (limit.nc) **Al1->restartdyn***
    920 ! --------------------------------------------------
    921 ! NB: List of the variables to be written in limit.nc
    922 !     (by writelim.F, subroutine of 1DUTILS.h):
    923 !        phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
    924 !        phy_fter,phy_foce,phy_flic,phy_fsic)
    925 !------------------------------------------------------------------------
    926       do i=1,yd
    927         phy_nat(i)  = nat_surf
    928         phy_alb(i)  = albedo
    929         phy_sst(i)  = tsurf ! read_tsurf1d will be used instead
    930         phy_rug(i)  = rugos
    931         phy_fter(i) = pctsrf(1,is_ter)
    932         phy_foce(i) = pctsrf(1,is_oce)
    933         phy_fsic(i) = pctsrf(1,is_sic)
    934         phy_flic(i) = pctsrf(1,is_lic)
    935       enddo
    936 
    937 ! fabrication de limit.nc
    938       call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,             &
    939      &               phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
    940 
    941 
    942       call phys_state_var_end
    943 !Al1
    944       if (restart) then
    945         print*,'call to restart dyn 1d'
    946         Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs,          &
    947      &              u,v,temp,q,omega2)
    948 
    949        print*,'fnday,annee_ref,day_ref,day_ini',                            &
    950      &     fnday,annee_ref,day_ref,day_ini
    951 !**      call ymds2ju(annee_ref,mois,day_ini,heure,day)
    952        day = day_ini
    953        day_end = day_ini + nday
    954        daytime = day + time_ini/24. ! 1st day and initial time of the simulation
    955 
    956 ! Print out the actual date of the beginning of the simulation :
    957        call ju2ymds(daytime, an, mois, jour, heure)
    958        print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
    959 
    960        day = int(daytime)
    961        time=daytime-day
    962  
    963        print*,'****** intialised fields from restart1dyn *******'
    964        print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
    965        print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
    966        print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
    967 ! raz for safety
    968        do l=1,llm
    969          dq_dyn(l,1) = 0.
    970        enddo
    971       endif
    972 !Al1 ================  end restart =================================
    973       IF (ecrit_slab_oc.eq.1) then
    974          open(97,file='div_slab.dat',STATUS='UNKNOWN')
    975        elseif (ecrit_slab_oc.eq.0) then
    976          open(97,file='div_slab.dat',STATUS='OLD')
    977        endif
    978 !
    979 !---------------------------------------------------------------------
    980 !    Initialize target profile for RHT nudging if needed
    981 !---------------------------------------------------------------------
    982       if (nudge(inudge_RHT)) then
    983         call nudge_RHT_init(plev,play,temp,q(:,1),t_targ,rh_targ)
    984       endif
    985       if (nudge(inudge_UV)) then
    986         call nudge_UV_init(plev,play,u,v,u_targ,v_targ)
    987       endif
    988 !
    989 !=====================================================================
    990        CALL iophys_ini
    991 ! START OF THE TEMPORAL LOOP :
    992 !=====================================================================
    993            
    994       it_end = nint(fnday*day_step)
    995 !test JLD     it_end = 10
    996       do while(it.le.it_end)
    997 
    998        if (prt_level.ge.1) then
    999          print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',                       &
    1000      &             it,day,time,it_end,day_step
    1001          print*,'PAS DE TEMPS ',timestep
    1002        endif
    1003 !Al1 demande de restartphy.nc
    1004        if (it.eq.it_end) lastcall=.True.
    1005 
    1006 !---------------------------------------------------------------------
    1007 ! Interpolation of forcings in time and onto model levels
    1008 !---------------------------------------------------------------------
    1009 
    1010 #include "1D_interp_cases.h"
    1011 
    1012       if (forcing_GCM2SCM) then
    1013         write (*,*) 'forcing_GCM2SCM not yet implemented'
    1014         stop 'in time loop'
    1015       endif ! forcing_GCM2SCM
    1016 
    1017 !---------------------------------------------------------------------
    1018 !  Geopotential :
    1019 !---------------------------------------------------------------------
    1020 
    1021         phi(1)=RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
    1022         do l = 1, llm-1
    1023           phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*                           &
    1024      &    (play(l)-play(l+1))/(play(l)+play(l+1))
    1025         enddo
    1026 
    1027 !---------------------------------------------------------------------
    1028 ! Listing output for debug prt_level>=1
    1029 !---------------------------------------------------------------------
    1030        if (prt_level>=1) then
    1031          print *,' avant physiq : -------- day time ',day,time
    1032          write(*,*) 'firstcall,lastcall,phis',                               &
    1033      &               firstcall,lastcall,phis
    1034        end if
    1035        if (prt_level>=5) then
    1036          write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l',                   &
    1037      &        'presniv','plev','play','phi'
    1038          write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l,                   &
    1039      &         presnivs(l),plev(l),play(l),phi(l),l=1,llm)
    1040          write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l',                    &
    1041      &         'presniv','u','v','temp','q1','q2','omega2'
    1042          write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l,         &
    1043      &   presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
    1044        endif
    1045 
    1046 !---------------------------------------------------------------------
    1047 !   Call physiq :
    1048 !---------------------------------------------------------------------
    1049        call physiq(ngrid,llm, &
    1050                     firstcall,lastcall,timestep, &
    1051                     plev,play,phi,phis,presnivs, &
    1052                     u,v, rot, temp,q,omega2, &
    1053                     du_phys,dv_phys,dt_phys,dq,dpsrf)
    1054                 firstcall=.false.
    1055 
    1056 !---------------------------------------------------------------------
    1057 ! Listing output for debug
    1058 !---------------------------------------------------------------------
    1059         if (prt_level>=5) then
    1060           write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l',                  &
    1061      &        'presniv','plev','play','phi'
    1062           write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l,                  &
    1063      &    presnivs(l),plev(l),play(l),phi(l),l=1,llm)
    1064           write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l',                   &
    1065      &         'presniv','u','v','temp','q1','q2','omega2'
    1066           write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l,       &
    1067      &    presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
    1068           write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l',                   &
    1069      &         'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'   
    1070            write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l,            &
    1071      &      presnivs(l),86400*du_phys(l),86400*dv_phys(l),                   &
    1072      &       86400*dt_phys(l),86400*dq(l,1),dq(l,2),l=1,llm)
    1073           write(*,*) 'dpsrf',dpsrf
    1074         endif
    1075 !---------------------------------------------------------------------
    1076 !   Add physical tendencies :
    1077 !---------------------------------------------------------------------
    1078 
    1079        fcoriolis=2.*sin(rpi*xlat/180.)*romega
    1080        if (forcing_radconv .or. forcing_fire) then
    1081          fcoriolis=0.0
    1082          dt_cooling=0.0
    1083          d_t_adv=0.0
    1084          d_q_adv=0.0
    1085        endif
    1086 !      print*, 'calcul de fcoriolis ', fcoriolis
    1087 
    1088        if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice            &
    1089      &    .or.forcing_amma .or. forcing_type.eq.101) then
    1090          fcoriolis=0.0 ; ug=0. ; vg=0.
    1091        endif
    1092 
    1093        if(forcing_rico) then
    1094           dt_cooling=0.
    1095        endif
    1096 
    1097 !CRio:Attention modif sp??cifique cas de Caroline
    1098       if (forcing_type==-1) then
    1099          fcoriolis=0.
    1100 !Nudging
    1101        
    1102 !on calcule dt_cooling
    1103         do l=1,llm
    1104         if (play(l).ge.20000.) then
    1105             dt_cooling(l)=-1.5/86400.
    1106         elseif ((play(l).ge.10000.).and.((play(l).lt.20000.))) then
    1107             dt_cooling(l)=-1.5/86400.*(play(l)-10000.)/(10000.)-1./86400.*(20000.-play(l))/10000.*(temp(l)-200.)
    1108         else
    1109             dt_cooling(l)=-1.*(temp(l)-200.)/86400.
    1110         endif
    1111         enddo
    1112 
    1113       endif     
    1114 !RC
    1115       if (forcing_sandu) then
    1116          ug(1:llm)=u_mod(1:llm)
    1117          vg(1:llm)=v_mod(1:llm)
    1118       endif
    1119 
    1120       IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
    1121                                    fcoriolis, xlat,mxcalc
    1122 
    1123 !       print *,'u-ug=',u-ug
    1124 
    1125 !!!!!!!!!!!!!!!!!!!!!!!!
    1126 ! Geostrophic wind
    1127 ! Le calcul ci dessous est insuffisamment precis
    1128 !      du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
    1129 !      dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
    1130 !!!!!!!!!!!!!!!!!!!!!!!!
    1131        sfdt = sin(0.5*fcoriolis*timestep)
    1132        cfdt = cos(0.5*fcoriolis*timestep)
    1133 !       print *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep
    1134 !
    1135         du_age(1:mxcalc)= -2.*sfdt/timestep*                                &
    1136      &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
    1137      &           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
    1138 !!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
    1139 !
    1140        dv_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
    1141      &          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +                           &
    1142      &           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
    1143 !!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
    1144 !
    1145 !!!!!!!!!!!!!!!!!!!!!!!!
    1146 !  Nudging
    1147 !!!!!!!!!!!!!!!!!!!!!!!!
    1148       d_t_nudge(:) = 0.
    1149       d_q_nudge(:,:) = 0.
    1150       d_u_nudge(:) = 0.
    1151       d_v_nudge(:) = 0.
    1152       if (nudge(inudge_RHT)) then
    1153         call nudge_RHT(timestep,plev,play,t_targ,rh_targ,temp,q(:,1),     &
    1154     &                  d_t_nudge,d_q_nudge(:,1))
    1155       endif
    1156       if (nudge(inudge_UV)) then
    1157         call nudge_UV(timestep,plev,play,u_targ,v_targ,u,v,     &
    1158     &                  d_u_nudge,d_v_nudge)
    1159       endif
    1160 !
    1161 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1162 !         call  writefield_phy('dv_age' ,dv_age,llm)
    1163 !         call  writefield_phy('du_age' ,du_age,llm)
    1164 !         call  writefield_phy('du_phys' ,du_phys,llm)
    1165 !         call  writefield_phy('u_tend' ,u,llm)
    1166 !         call  writefield_phy('u_g' ,ug,llm)
    1167 !
    1168 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1169 !! Increment state variables
    1170 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1171     IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
    1172 
    1173 ! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h
    1174 ! au dessus de 700hpa, on relaxe vers les profils initiaux
    1175       if (forcing_sandu .OR. forcing_astex) then
    1176 #include "1D_nudge_sandu_astex.h"
    1177       else
    1178         u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
    1179      &              du_phys(1:mxcalc)                                       &
    1180      &             +du_age(1:mxcalc)+du_adv(1:mxcalc)                       &
    1181      &             +d_u_nudge(1:mxcalc) )           
    1182         v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
    1183      &              dv_phys(1:mxcalc)                                       &
    1184      &             +dv_age(1:mxcalc)+dv_adv(1:mxcalc)                       &
    1185      &             +d_v_nudge(1:mxcalc) )
    1186         q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
    1187      &                dq(1:mxcalc,:)                                        &
    1188      &               +d_q_adv(1:mxcalc,:)                                   &
    1189      &               +d_q_nudge(1:mxcalc,:) )
    1190 
    1191         if (prt_level.ge.3) then
    1192           print *,                                                          &
    1193      &    'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ',         &
    1194      &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
    1195            print* ,'dv_phys=',dv_phys
    1196            print* ,'dv_age=',dv_age
    1197            print* ,'dv_adv=',dv_adv
    1198            print* ,'d_v_nudge=',d_v_nudge
    1199            print*, v
    1200            print*, vg
    1201         endif
    1202 
    1203         temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
    1204      &              dt_phys(1:mxcalc)                                       &
    1205      &             +d_t_adv(1:mxcalc)                                      &
    1206      &             +d_t_nudge(1:mxcalc)                                      &
    1207      &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
    1208 
    1209       endif  ! forcing_sandu or forcing_astex
    1210 
    1211         teta=temp*(pzero/play)**rkappa
    1212 !
    1213 !---------------------------------------------------------------------
    1214 !   Nudge soil temperature if requested
    1215 !---------------------------------------------------------------------
    1216 
    1217       IF (nudge_tsoil .AND. .NOT. lastcall) THEN
    1218        ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:)                     &
    1219      &  -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
    1220       ENDIF
    1221 
    1222 !---------------------------------------------------------------------
    1223 !   Add large-scale tendencies (advection, etc) :
    1224 !---------------------------------------------------------------------
    1225 
    1226 !cc nrlmd
    1227 !cc        tmpvar=teta
    1228 !cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
    1229 !cc
    1230 !cc        teta(1:mxcalc)=tmpvar(1:mxcalc)
    1231 !cc        tmpvar(:)=q(:,1)
    1232 !cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
    1233 !cc        q(1:mxcalc,1)=tmpvar(1:mxcalc)
    1234 !cc        tmpvar(:)=q(:,2)
    1235 !cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
    1236 !cc        q(1:mxcalc,2)=tmpvar(1:mxcalc)
    1237 
    1238    END IF ! end if tendency of tendency should be added
    1239 
    1240 !---------------------------------------------------------------------
    1241 !   Air temperature :
    1242 !---------------------------------------------------------------------       
    1243         if (lastcall) then
    1244           print*,'Pas de temps final ',it
    1245           call ju2ymds(daytime, an, mois, jour, heure)
    1246           print*,'a la date : a m j h',an, mois, jour ,heure/3600.
    1247         endif
    1248 
    1249 !  incremente day time
    1250 !        print*,'daytime bef',daytime,1./day_step
    1251         daytime = daytime+1./day_step
    1252 !Al1dbg
    1253         day = int(daytime+0.1/day_step)
    1254 !        time = max(daytime-day,0.0)
    1255 !Al1&jyg: correction de bug
    1256 !cc        time = real(mod(it,day_step))/day_step
    1257         time = time_ini/24.+real(mod(it,day_step))/day_step
    1258 !        print*,'daytime nxt time',daytime,time
    1259         it=it+1
    1260 
    1261       enddo
    1262 
    1263 !Al1
    1264       if (ecrit_slab_oc.ne.-1) close(97)
    1265 
    1266 !Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
    1267 ! -------------------------------------
    1268        call dyn1dredem("restart1dyn.nc",                                    &
    1269      &              plev,play,phi,phis,presnivs,                            &
    1270      &              u,v,temp,q,omega2)
    1271 
    1272         CALL abort_gcm ('lmdz1d   ','The End  ',0)
    1273 
    1274       end
    127527
    127628#include "1DUTILS.h"
Note: See TracChangeset for help on using the changeset viewer.