Changeset 2126 for LMDZ5/trunk/libf


Ignore:
Timestamp:
Oct 15, 2014, 2:03:57 AM (10 years ago)
Author:
fhourdin
Message:

Introduction du cas Dice couplé atmosphere/surface
+ nouveau paramètre de contrôle f_ri_cd_min, seuil minimum
sur la fonction f(Ri) en facteur du coefficient de traîné neutre.
Par défaut : f_ri_cd_min=0.1 (comme avant)

Introduction of the coupled atmosphere/surface dice case.
+ new parameter f_ri_cd_min, minimum threshold on the f(Ri) fonction
that multiply the neutral drag coefficient.

Location:
LMDZ5/trunk/libf/phylmd
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/1DUTILS.h

    r2117 r2126  
    140140       CALL getin('ok_flux_surf',ok_flux_surf)
    141141
     142!Config  Key  = ok_prescr_ust
     143!Config  Desc = ustar impose ou non
     144!Config  Def  = false
     145!Config  Help = ustar impose ou non
     146       ok_prescr_ust = .false.
     147       CALL getin('ok_prescr_ust',ok_prescr_ust)
     148
    142149!Config  Key  = ok_old_disvert
    143150!Config  Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
    … …  
    24132420          end
    24142421 
     2422!=====================================================================
     2423       SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof   &
     2424     &         ,th_prof,qv_prof,u_prof,v_prof,o3_prof                     &
     2425     &         ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof         &
     2426     &         ,th_mod,qv_mod,u_mod,v_mod,o3_mod                          &
     2427     &         ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
     2428 
     2429       implicit none
     2430 
     2431#include "dimensions.h"
     2432
     2433!-------------------------------------------------------------------------
     2434! Vertical interpolation of Dice forcing data onto model levels
     2435!-------------------------------------------------------------------------
     2436 
     2437       integer nlevmax
     2438       parameter (nlevmax=41)
     2439       integer nlev_dice,mxcalc,nt_dice
     2440 
     2441       real play(llm), plev_prof(nlev_dice)
     2442       real th_prof(nlev_dice),qv_prof(nlev_dice)
     2443       real u_prof(nlev_dice),v_prof(nlev_dice)
     2444       real o3_prof(nlev_dice)
     2445       real ht_prof(nlev_dice),hq_prof(nlev_dice)
     2446       real hu_prof(nlev_dice),hv_prof(nlev_dice)
     2447       real w_prof(nlev_dice),omega_prof(nlev_dice)
     2448 
     2449       real th_mod(llm),qv_mod(llm)
     2450       real u_mod(llm),v_mod(llm), o3_mod(llm)
     2451       real ht_mod(llm),hq_mod(llm)
     2452       real hu_mod(llm),hv_mod(llm),w_mod(llm),omega_mod(llm)
     2453 
     2454       integer l,k,k1,k2,kp
     2455       real aa,frac,frac1,frac2,fact
     2456 
     2457       do l = 1, llm
     2458
     2459        if (play(l).ge.plev_prof(nlev_dice)) then
     2460 
     2461        mxcalc=l
     2462         k1=0
     2463         k2=0
     2464
     2465         if (play(l).le.plev_prof(1)) then
     2466
     2467         do k = 1, nlev_dice-1
     2468          if (play(l).le.plev_prof(k) .and. play(l).gt.plev_prof(k+1)) then
     2469            k1=k
     2470            k2=k+1
     2471          endif
     2472         enddo
     2473
     2474         if (k1.eq.0 .or. k2.eq.0) then
     2475          write(*,*) 'PB! k1, k2 = ',k1,k2
     2476          write(*,*) 'l,play(l) = ',l,play(l)/100
     2477         do k = 1, nlev_dice-1
     2478          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
     2479         enddo
     2480         endif
     2481
     2482         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
     2483         th_mod(l)= th_prof(k2) - frac*(th_prof(k2)-th_prof(k1))
     2484         qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
     2485         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
     2486         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
     2487         o3_mod(l)= o3_prof(k2) - frac*(o3_prof(k2)-o3_prof(k1))
     2488         ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
     2489         hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
     2490         hu_mod(l)= hu_prof(k2) - frac*(hu_prof(k2)-hu_prof(k1))
     2491         hv_mod(l)= hv_prof(k2) - frac*(hv_prof(k2)-hv_prof(k1))
     2492         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
     2493         omega_mod(l)= omega_prof(k2) - frac*(omega_prof(k2)-omega_prof(k1))
     2494     
     2495         else !play>plev_prof(1)
     2496
     2497         k1=1
     2498         k2=2
     2499         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
     2500         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
     2501         th_mod(l)= frac1*th_prof(k1) - frac2*th_prof(k2)
     2502         qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2)
     2503         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
     2504         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
     2505         o3_mod(l)= frac1*o3_prof(k1) - frac2*o3_prof(k2)
     2506         ht_mod(l)= frac1*ht_prof(k1) - frac2*ht_prof(k2)
     2507         hq_mod(l)= frac1*hq_prof(k1) - frac2*hq_prof(k2)
     2508         hu_mod(l)= frac1*hu_prof(k1) - frac2*hu_prof(k2)
     2509         hv_mod(l)= frac1*hv_prof(k1) - frac2*hv_prof(k2)
     2510         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
     2511         omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2)
     2512
     2513         endif ! play.le.plev_prof(1)
     2514
     2515        else ! above max altitude of forcing file
     2516 
     2517!jyg
     2518         fact=20.*(plev_prof(nlev_dice)-play(l))/plev_prof(nlev_dice) !jyg
     2519         fact = max(fact,0.)                                           !jyg
     2520         fact = exp(-fact)                                             !jyg
     2521         th_mod(l)= th_prof(nlev_dice)                                 !jyg
     2522         qv_mod(l)= qv_prof(nlev_dice)*fact                            !jyg
     2523         u_mod(l)= u_prof(nlev_dice)*fact                              !jyg
     2524         v_mod(l)= v_prof(nlev_dice)*fact                              !jyg
     2525         o3_mod(l)= o3_prof(nlev_dice)*fact                            !jyg
     2526         ht_mod(l)= ht_prof(nlev_dice)                                 !jyg
     2527         hq_mod(l)= hq_prof(nlev_dice)*fact                            !jyg
     2528         hu_mod(l)= hu_prof(nlev_dice)                                 !jyg
     2529         hv_mod(l)= hv_prof(nlev_dice)                                 !jyg
     2530         w_mod(l)= 0.                                                  !jyg
     2531         omega_mod(l)= 0.                                              !jyg
     2532 
     2533        endif ! play
     2534 
     2535       enddo ! l
     2536
     2537!       do l = 1,llm
     2538!       print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
     2539!     $        l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
     2540!       enddo
     2541 
     2542          return
     2543          end
     2544
    24152545!======================================================================
    24162546        SUBROUTINE interp_astex_time(day,day1,annee_ref                    &
    … …  
    26452775
    26462776!======================================================================
     2777        SUBROUTINE interp_dice_time(day,day1,annee_ref                    &
     2778     &             ,year_ini_dice,day_ini_dice,nt_dice,dt_dice            &
     2779     &             ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice       &
     2780     &             ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice         &
     2781     &             ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice     &
     2782     &             ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof         &
     2783     &             ,ustar_prof,psurf_prof,ug_prof,vg_prof                 &
     2784     &             ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof)
     2785        implicit none
     2786
     2787!---------------------------------------------------------------------------------------
     2788! Time interpolation of a 2D field to the timestep corresponding to day
     2789!
     2790! day: current julian day (e.g. 717538.2)
     2791! day1: first day of the simulation
     2792! nt_dice: total nb of data in the forcing (e.g. 145 for Dice)
     2793! dt_dice: total time interval (in sec) between 2 forcing data (e.g. 30min. for Dice)
     2794!---------------------------------------------------------------------------------------
     2795
     2796#include "compar1d.h"
     2797
     2798! inputs:
     2799        integer annee_ref
     2800        integer nt_dice,nlev_dice
     2801        integer year_ini_dice
     2802        real day, day1,day_ini_dice,dt_dice
     2803        real shf_dice(nt_dice),lhf_dice(nt_dice),lwup_dice(nt_dice)
     2804        real swup_dice(nt_dice),tg_dice(nt_dice),ustar_dice(nt_dice)
     2805        real psurf_dice(nt_dice),ug_dice(nt_dice),vg_dice(nt_dice)
     2806        real ht_dice(nlev_dice,nt_dice),hq_dice(nlev_dice,nt_dice)
     2807        real hu_dice(nlev_dice,nt_dice),hv_dice(nlev_dice,nt_dice)
     2808        real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice)
     2809! outputs:
     2810        real tg_prof,shf_prof,lhf_prof,lwup_prof,swup_prof
     2811        real ustar_prof,psurf_prof,ug_prof,vg_prof
     2812        real ht_prof(nlev_dice),hq_prof(nlev_dice)
     2813        real hu_prof(nlev_dice),hv_prof(nlev_dice)
     2814        real w_prof(nlev_dice),omega_prof(nlev_dice)
     2815! local:
     2816        integer it_dice1, it_dice2,k
     2817        real timeit,time_dice1,time_dice2,frac
     2818
     2819
     2820        if (forcing_type.eq.7) then
     2821! Check that initial day of the simulation consistent with Dice period:
     2822       print *,'annee_ref=',annee_ref
     2823       print *,'day1=',day1
     2824       print *,'day_ini_dice=',day_ini_dice
     2825       if (annee_ref.ne.1999) then
     2826        print*,'Pour Dice, annee_ref doit etre 1999'
     2827        print*,'Changer annee_ref dans run.def'
     2828        stop
     2829       endif
     2830       if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice) then
     2831        print*,'Dice a debute le 23 Oct 1999 (jour julien=296)'
     2832        print*,'Changer dayref dans run.def',day1,day_ini_dice
     2833        stop
     2834       endif
     2835       if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice+2) then
     2836        print*,'Dice a fini le 25 Oct 1999 (jour julien=298)'
     2837        print*,'Changer dayref ou nday dans run.def',day1,day_ini_dice
     2838        stop
     2839       endif
     2840
     2841       endif
     2842
     2843! Determine timestep relative to the 1st day of TOGA-COARE:
     2844!       timeit=(day-day1)*86400.
     2845!       if (annee_ref.eq.1992) then
     2846!        timeit=(day-day_ini_dice)*86400.
     2847!       else
     2848!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
     2849!       endif
     2850      timeit=(day-day_ini_dice)*86400
     2851
     2852! Determine the closest observation times:
     2853       it_dice1=INT(timeit/dt_dice)+1
     2854       it_dice2=it_dice1 + 1
     2855       time_dice1=(it_dice1-1)*dt_dice
     2856       time_dice2=(it_dice2-1)*dt_dice
     2857
     2858       if (it_dice1 .ge. nt_dice) then
     2859        write(*,*) 'PB-stop: day, it_dice1, it_dice2, timeit: ',day,it_dice1,it_dice2,timeit/86400.
     2860        stop
     2861       endif
     2862
     2863! time interpolation:
     2864       frac=(time_dice2-timeit)/(time_dice2-time_dice1)
     2865       frac=max(frac,0.0)
     2866
     2867       shf_prof = shf_dice(it_dice2)-frac*(shf_dice(it_dice2)-shf_dice(it_dice1))
     2868       lhf_prof = lhf_dice(it_dice2)-frac*(lhf_dice(it_dice2)-lhf_dice(it_dice1))
     2869       lwup_prof = lwup_dice(it_dice2)-frac*(lwup_dice(it_dice2)-lwup_dice(it_dice1))
     2870       swup_prof = swup_dice(it_dice2)-frac*(swup_dice(it_dice2)-swup_dice(it_dice1))
     2871       tg_prof = tg_dice(it_dice2)-frac*(tg_dice(it_dice2)-tg_dice(it_dice1))
     2872       ustar_prof = ustar_dice(it_dice2)-frac*(ustar_dice(it_dice2)-ustar_dice(it_dice1))
     2873       psurf_prof = psurf_dice(it_dice2)-frac*(psurf_dice(it_dice2)-psurf_dice(it_dice1))
     2874       ug_prof = ug_dice(it_dice2)-frac*(ug_dice(it_dice2)-ug_dice(it_dice1))
     2875       vg_prof = vg_dice(it_dice2)-frac*(vg_dice(it_dice2)-vg_dice(it_dice1))
     2876
     2877!        print*,
     2878!     :'day,annee_ref,day_ini_dice,timeit,it_dice1,it_dice2,SST:',
     2879!     :day,annee_ref,day_ini_dice,timeit/86400.,it_dice1,it_dice2,ts_prof
     2880
     2881       do k=1,nlev_dice
     2882        ht_prof(k) = ht_dice(k,it_dice2)-frac*(ht_dice(k,it_dice2)-ht_dice(k,it_dice1))
     2883        hq_prof(k) = hq_dice(k,it_dice2)-frac*(hq_dice(k,it_dice2)-hq_dice(k,it_dice1))
     2884        hu_prof(k) = hu_dice(k,it_dice2)-frac*(hu_dice(k,it_dice2)-hu_dice(k,it_dice1))
     2885        hv_prof(k) = hv_dice(k,it_dice2)-frac*(hv_dice(k,it_dice2)-hv_dice(k,it_dice1))
     2886        w_prof(k) = w_dice(k,it_dice2)-frac*(w_dice(k,it_dice2)-w_dice(k,it_dice1))
     2887        omega_prof(k) = omega_dice(k,it_dice2)-frac*(omega_dice(k,it_dice2)-omega_dice(k,it_dice1))
     2888        enddo
     2889
     2890        return
     2891        END
     2892
     2893!======================================================================
    26472894        SUBROUTINE interp_armcu_time(day,day1,annee_ref                    &
    26482895     &             ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu         &
    … …  
    31723419         return
    31733420         end subroutine read_fire
     3421!=====================================================================
     3422      subroutine read_dice(fich_dice,nlevel,ntime                         &
     3423     &     ,zz,pres,th,qv,u,v,o3                                          &
     3424     &     ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg                        &
     3425     &     ,hadvt,hadvq,hadvu,hadvv,w,omega)
     3426
     3427!program reading initial profils and forcings of the Dice case study
     3428
     3429
     3430      implicit none
     3431
     3432#include "netcdf.inc"
     3433
     3434      integer ntime,nlevel
     3435      integer l,k
     3436      character*80 :: fich_dice
     3437      real*8 time(ntime)
     3438      real*8 zz(nlevel)
     3439
     3440      real*8 th(nlevel),pres(nlevel)
     3441      real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel)
     3442      real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime)
     3443      real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime)
     3444      real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime)
     3445      real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime)
     3446
     3447      integer nid, ierr
     3448      integer nbvar3d
     3449      parameter(nbvar3d=30)
     3450      integer var3didin(nbvar3d)
     3451
     3452      ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid)
     3453      if (ierr.NE.NF_NOERR) then
     3454         write(*,*) 'ERROR: Pb opening forcings nc file '
     3455         write(*,*) NF_STRERROR(ierr)
     3456         stop ""
     3457      endif
     3458
     3459
     3460       ierr=NF_INQ_VARID(nid,"height",var3didin(1))
     3461         if(ierr/=NF_NOERR) then
     3462           write(*,*) NF_STRERROR(ierr)
     3463           stop 'height'
     3464         endif
     3465
     3466       ierr=NF_INQ_VARID(nid,"pf",var3didin(11))
     3467         if(ierr/=NF_NOERR) then
     3468           write(*,*) NF_STRERROR(ierr)
     3469           stop 'pf'
     3470         endif
     3471
     3472      ierr=NF_INQ_VARID(nid,"theta",var3didin(12))
     3473         if(ierr/=NF_NOERR) then
     3474           write(*,*) NF_STRERROR(ierr)
     3475           stop 'theta'
     3476         endif
     3477
     3478      ierr=NF_INQ_VARID(nid,"qv",var3didin(13))
     3479         if(ierr/=NF_NOERR) then
     3480           write(*,*) NF_STRERROR(ierr)
     3481           stop 'qv'
     3482         endif
     3483
     3484      ierr=NF_INQ_VARID(nid,"u",var3didin(14))
     3485         if(ierr/=NF_NOERR) then
     3486           write(*,*) NF_STRERROR(ierr)
     3487           stop 'u'
     3488         endif
     3489
     3490      ierr=NF_INQ_VARID(nid,"v",var3didin(15))
     3491         if(ierr/=NF_NOERR) then
     3492           write(*,*) NF_STRERROR(ierr)
     3493           stop 'v'
     3494         endif
     3495
     3496      ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16))
     3497         if(ierr/=NF_NOERR) then
     3498           write(*,*) NF_STRERROR(ierr)
     3499           stop 'o3'
     3500         endif
     3501
     3502      ierr=NF_INQ_VARID(nid,"shf",var3didin(2))
     3503         if(ierr/=NF_NOERR) then
     3504           write(*,*) NF_STRERROR(ierr)
     3505           stop 'shf'
     3506         endif
     3507
     3508      ierr=NF_INQ_VARID(nid,"lhf",var3didin(3))
     3509         if(ierr/=NF_NOERR) then
     3510           write(*,*) NF_STRERROR(ierr)
     3511           stop 'lhf'
     3512         endif
     3513     
     3514      ierr=NF_INQ_VARID(nid,"lwup",var3didin(4))
     3515         if(ierr/=NF_NOERR) then
     3516           write(*,*) NF_STRERROR(ierr)
     3517           stop 'lwup'
     3518         endif
     3519
     3520      ierr=NF_INQ_VARID(nid,"swup",var3didin(5))
     3521         if(ierr/=NF_NOERR) then
     3522           write(*,*) NF_STRERROR(ierr)
     3523           stop 'dqtdx'
     3524         endif
     3525
     3526      ierr=NF_INQ_VARID(nid,"Tg",var3didin(6))
     3527         if(ierr/=NF_NOERR) then
     3528           write(*,*) NF_STRERROR(ierr)
     3529           stop 'Tg'
     3530      endif
     3531
     3532      ierr=NF_INQ_VARID(nid,"ustar",var3didin(7))
     3533         if(ierr/=NF_NOERR) then
     3534           write(*,*) NF_STRERROR(ierr)
     3535           stop 'ustar'
     3536      endif
     3537
     3538      ierr=NF_INQ_VARID(nid,"psurf",var3didin(8))
     3539         if(ierr/=NF_NOERR) then
     3540           write(*,*) NF_STRERROR(ierr)
     3541           stop 'psurf'
     3542      endif
     3543
     3544      ierr=NF_INQ_VARID(nid,"Ug",var3didin(9))
     3545         if(ierr/=NF_NOERR) then
     3546           write(*,*) NF_STRERROR(ierr)
     3547           stop 'Ug'
     3548      endif
     3549
     3550      ierr=NF_INQ_VARID(nid,"Vg",var3didin(10))
     3551         if(ierr/=NF_NOERR) then
     3552           write(*,*) NF_STRERROR(ierr)
     3553           stop 'Vg'
     3554      endif
     3555
     3556      ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17))
     3557         if(ierr/=NF_NOERR) then
     3558           write(*,*) NF_STRERROR(ierr)
     3559           stop 'hadvT'
     3560      endif
     3561
     3562      ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18))
     3563         if(ierr/=NF_NOERR) then
     3564           write(*,*) NF_STRERROR(ierr)
     3565           stop 'hadvq'
     3566      endif
     3567
     3568      ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19))
     3569         if(ierr/=NF_NOERR) then
     3570           write(*,*) NF_STRERROR(ierr)
     3571           stop 'hadvu'
     3572      endif
     3573
     3574      ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20))
     3575         if(ierr/=NF_NOERR) then
     3576           write(*,*) NF_STRERROR(ierr)
     3577           stop 'hadvv'
     3578      endif
     3579
     3580      ierr=NF_INQ_VARID(nid,"w",var3didin(21))
     3581         if(ierr/=NF_NOERR) then
     3582           write(*,*) NF_STRERROR(ierr)
     3583           stop 'w'
     3584      endif
     3585
     3586      ierr=NF_INQ_VARID(nid,"omega",var3didin(22))
     3587         if(ierr/=NF_NOERR) then
     3588           write(*,*) NF_STRERROR(ierr)
     3589           stop 'omega'
     3590      endif
     3591!dimensions lecture
     3592!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
     3593 
     3594#ifdef NC_DOUBLE
     3595         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
     3596#else
     3597         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
     3598#endif
     3599         if(ierr/=NF_NOERR) then
     3600            write(*,*) NF_STRERROR(ierr)
     3601            stop "getvarup"
     3602         endif
     3603!          write(*,*)'lecture zz ok',zz
     3604 
     3605#ifdef NC_DOUBLE
     3606         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pres)
     3607#else
     3608         ierr = NF_GET_VAR_REAL(nid,var3didin(11),pres)
     3609#endif
     3610         if(ierr/=NF_NOERR) then
     3611            write(*,*) NF_STRERROR(ierr)
     3612            stop "getvarup"
     3613         endif
     3614!          write(*,*)'lecture pres ok',pres
     3615
     3616#ifdef NC_DOUBLE
     3617         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),th)
     3618#else
     3619         ierr = NF_GET_VAR_REAL(nid,var3didin(12),th)
     3620#endif
     3621         if(ierr/=NF_NOERR) then
     3622            write(*,*) NF_STRERROR(ierr)
     3623            stop "getvarup"
     3624         endif
     3625!          write(*,*)'lecture th ok',th
     3626
     3627#ifdef NC_DOUBLE
     3628         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),qv)
     3629#else
     3630         ierr = NF_GET_VAR_REAL(nid,var3didin(13),qv)
     3631#endif
     3632         if(ierr/=NF_NOERR) then
     3633            write(*,*) NF_STRERROR(ierr)
     3634            stop "getvarup"
     3635         endif
     3636!          write(*,*)'lecture qv ok',qv
     3637 
     3638#ifdef NC_DOUBLE
     3639         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),u)
     3640#else
     3641         ierr = NF_GET_VAR_REAL(nid,var3didin(14),u)
     3642#endif
     3643         if(ierr/=NF_NOERR) then
     3644            write(*,*) NF_STRERROR(ierr)
     3645            stop "getvarup"
     3646         endif
     3647!          write(*,*)'lecture u ok',u
     3648
     3649#ifdef NC_DOUBLE
     3650         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),v)
     3651#else
     3652         ierr = NF_GET_VAR_REAL(nid,var3didin(15),v)
     3653#endif
     3654         if(ierr/=NF_NOERR) then
     3655            write(*,*) NF_STRERROR(ierr)
     3656            stop "getvarup"
     3657         endif
     3658!          write(*,*)'lecture v ok',v
     3659
     3660#ifdef NC_DOUBLE
     3661         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),o3)
     3662#else
     3663         ierr = NF_GET_VAR_REAL(nid,var3didin(16),o3)
     3664#endif
     3665         if(ierr/=NF_NOERR) then
     3666            write(*,*) NF_STRERROR(ierr)
     3667            stop "getvarup"
     3668         endif
     3669!          write(*,*)'lecture o3 ok',o3
     3670
     3671#ifdef NC_DOUBLE
     3672         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),shf)
     3673#else
     3674         ierr = NF_GET_VAR_REAL(nid,var3didin(2),shf)
     3675#endif
     3676         if(ierr/=NF_NOERR) then
     3677            write(*,*) NF_STRERROR(ierr)
     3678            stop "getvarup"
     3679         endif
     3680!          write(*,*)'lecture shf ok',shf
     3681
     3682#ifdef NC_DOUBLE
     3683         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),lhf)
     3684#else
     3685         ierr = NF_GET_VAR_REAL(nid,var3didin(3),lhf)
     3686#endif
     3687         if(ierr/=NF_NOERR) then
     3688            write(*,*) NF_STRERROR(ierr)
     3689            stop "getvarup"
     3690         endif
     3691!          write(*,*)'lecture lhf ok',lhf
     3692
     3693#ifdef NC_DOUBLE
     3694         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),lwup)
     3695#else
     3696         ierr = NF_GET_VAR_REAL(nid,var3didin(4),lwup)
     3697#endif
     3698         if(ierr/=NF_NOERR) then
     3699            write(*,*) NF_STRERROR(ierr)
     3700            stop "getvarup"
     3701         endif
     3702!          write(*,*)'lecture lwup ok',lwup
     3703
     3704#ifdef NC_DOUBLE
     3705         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),swup)
     3706#else
     3707         ierr = NF_GET_VAR_REAL(nid,var3didin(5),swup)
     3708#endif
     3709         if(ierr/=NF_NOERR) then
     3710            write(*,*) NF_STRERROR(ierr)
     3711            stop "getvarup"
     3712         endif
     3713!          write(*,*)'lecture swup ok',swup
     3714
     3715#ifdef NC_DOUBLE
     3716         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tg)
     3717#else
     3718         ierr = NF_GET_VAR_REAL(nid,var3didin(6),tg)
     3719#endif
     3720         if(ierr/=NF_NOERR) then
     3721            write(*,*) NF_STRERROR(ierr)
     3722            stop "getvarup"
     3723         endif
     3724!          write(*,*)'lecture tg ok',tg
     3725
     3726#ifdef NC_DOUBLE
     3727         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ustar)
     3728#else
     3729         ierr = NF_GET_VAR_REAL(nid,var3didin(7),ustar)
     3730#endif
     3731         if(ierr/=NF_NOERR) then
     3732            write(*,*) NF_STRERROR(ierr)
     3733            stop "getvarup"
     3734         endif
     3735!          write(*,*)'lecture ustar ok',ustar
     3736
     3737#ifdef NC_DOUBLE
     3738         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),psurf)
     3739#else
     3740         ierr = NF_GET_VAR_REAL(nid,var3didin(8),psurf)
     3741#endif
     3742         if(ierr/=NF_NOERR) then
     3743            write(*,*) NF_STRERROR(ierr)
     3744            stop "getvarup"
     3745         endif
     3746!          write(*,*)'lecture psurf ok',psurf
     3747
     3748#ifdef NC_DOUBLE
     3749         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),ug)
     3750#else
     3751         ierr = NF_GET_VAR_REAL(nid,var3didin(9),ug)
     3752#endif
     3753         if(ierr/=NF_NOERR) then
     3754            write(*,*) NF_STRERROR(ierr)
     3755            stop "getvarup"
     3756         endif
     3757!          write(*,*)'lecture ug ok',ug
     3758
     3759#ifdef NC_DOUBLE
     3760         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),vg)
     3761#else
     3762         ierr = NF_GET_VAR_REAL(nid,var3didin(10),vg)
     3763#endif
     3764         if(ierr/=NF_NOERR) then
     3765            write(*,*) NF_STRERROR(ierr)
     3766            stop "getvarup"
     3767         endif
     3768!          write(*,*)'lecture vg ok',vg
     3769
     3770#ifdef NC_DOUBLE
     3771         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hadvt)
     3772#else
     3773         ierr = NF_GET_VAR_REAL(nid,var3didin(17),hadvt)
     3774#endif
     3775         if(ierr/=NF_NOERR) then
     3776            write(*,*) NF_STRERROR(ierr)
     3777            stop "getvarup"
     3778         endif
     3779!          write(*,*)'lecture hadvt ok',hadvt
     3780
     3781#ifdef NC_DOUBLE
     3782         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),hadvq)
     3783#else
     3784         ierr = NF_GET_VAR_REAL(nid,var3didin(18),hadvq)
     3785#endif
     3786         if(ierr/=NF_NOERR) then
     3787            write(*,*) NF_STRERROR(ierr)
     3788            stop "getvarup"
     3789         endif
     3790!          write(*,*)'lecture hadvq ok',hadvq
     3791
     3792#ifdef NC_DOUBLE
     3793         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),hadvu)
     3794#else
     3795         ierr = NF_GET_VAR_REAL(nid,var3didin(19),hadvu)
     3796#endif
     3797         if(ierr/=NF_NOERR) then
     3798            write(*,*) NF_STRERROR(ierr)
     3799            stop "getvarup"
     3800         endif
     3801!          write(*,*)'lecture hadvu ok',hadvu
     3802
     3803#ifdef NC_DOUBLE
     3804         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),hadvv)
     3805#else
     3806         ierr = NF_GET_VAR_REAL(nid,var3didin(20),hadvv)
     3807#endif
     3808         if(ierr/=NF_NOERR) then
     3809            write(*,*) NF_STRERROR(ierr)
     3810            stop "getvarup"
     3811         endif
     3812!          write(*,*)'lecture hadvv ok',hadvv
     3813
     3814#ifdef NC_DOUBLE
     3815         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),w)
     3816#else
     3817         ierr = NF_GET_VAR_REAL(nid,var3didin(21),w)
     3818#endif
     3819         if(ierr/=NF_NOERR) then
     3820            write(*,*) NF_STRERROR(ierr)
     3821            stop "getvarup"
     3822         endif
     3823!          write(*,*)'lecture w ok',w
     3824
     3825#ifdef NC_DOUBLE
     3826         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),omega)
     3827#else
     3828         ierr = NF_GET_VAR_REAL(nid,var3didin(22),omega)
     3829#endif
     3830         if(ierr/=NF_NOERR) then
     3831            write(*,*) NF_STRERROR(ierr)
     3832            stop "getvarup"
     3833         endif
     3834!          write(*,*)'lecture omega ok',omega
     3835
     3836         return
     3837         end subroutine read_dice
    31743838!=====================================================================
    31753839!     Reads CIRC input files     

  • LMDZ5/trunk/libf/phylmd/1D_decl_cases.h

    r2117 r2126  
    3535        real u_mod(llm),v_mod(llm), ht_mod(llm),vt_mod(llm)
    3636        real hq_mod(llm),vq_mod(llm),qv_mod(llm),ql_mod(llm),qt_mod(llm)
     37        real th_mod(llm)
    3738
    3839        real ts_cur
    … …  
    9293        parameter (day_ini_fire=14)  ! 14 = 14Juil1987
    9394        parameter (heure_ini_fire=0.) !0h en secondes
     95
     96!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     97!Declarations specifiques au cas DICE     (MPL 02072013)
     98        character*80 :: fich_dice
     99        integer nlev_dice, nt_dice
     100        parameter (nlev_dice=70, nt_dice=145) 
     101        integer year_ini_dice, day_ini_dice, mth_ini_dice
     102        real heure_ini_dice
     103        real day_ju_ini_dice   ! Julian day of dice first day
     104        parameter (year_ini_dice=1999)
     105        parameter (mth_ini_dice=10)
     106        parameter (day_ini_dice=23)  ! 23 = 23 october 1999
     107        parameter (heure_ini_dice=68400.) !19UTC en secondes
     108        real dt_dice
     109        parameter (dt_dice=0.5*3600.) ! 1 forcage ttes les demi-heures
     110
     111!profils initiaux:
     112        real plev_dice(nlev_dice)
     113       
     114        real zz_dice(nlev_dice)
     115        real th_dice(nlev_dice),qv_dice(nlev_dice)
     116        real u_dice(nlev_dice), v_dice(nlev_dice),o3_dice(nlev_dice)
     117        real ht_dice(nlev_dice,nt_dice)
     118        real hq_dice(nlev_dice,nt_dice), hu_dice(nlev_dice,nt_dice)
     119        real hv_dice(nlev_dice,nt_dice)
     120        real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice)
     121        real o3_mod(llm),hu_mod(llm),hv_mod(llm)
     122        real th_dicei(nlev_dice),qv_dicei(nlev_dice)
     123        real u_dicei(nlev_dice), v_dicei(nlev_dice),o3_dicei(nlev_dice)
     124        real ht_dicei(nlev_dice)
     125        real hq_dicei(nlev_dice), hu_dicei(nlev_dice)
     126        real hv_dicei(nlev_dice)
     127        real w_dicei(nlev_dice),omega_dicei(nlev_dice)
     128
     129       
     130!forcings
     131        real shf_dice(nt_dice),lhf_dice(nt_dice)
     132        real lwup_dice(nt_dice),swup_dice(nt_dice)
     133        real tg_dice(nt_dice),ustar_dice(nt_dice),psurf_dice(nt_dice)
     134        real ug_dice(nt_dice),vg_dice(nt_dice)
     135
     136        real shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof
     137        real ustar_prof,psurf_prof,cdrag
     138        real ht_profd(nlev_dice),hq_profd(nlev_dice),hu_profd(nlev_dice)
     139        real hv_profd(nlev_dice),w_profd(nlev_dice)
     140        real omega_profd(nlev_dice),ug_profd,vg_profd
    94141
    95142!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    165212        real d_t_z(llm), d_q_z(llm)
    166213        real d_t_dyn_z(llm), d_q_dyn_z(llm)
     214        real d_u_z(llm),d_v_z(llm)
     215        real d_u_dyn(llm),d_v_dyn(llm)
     216        real d_u_dyn_z(llm),d_v_dyn_z(llm)
     217        real d_u_adv(llm),d_v_adv(llm)
    167218        real zz(llm)
    168219        real zfact

  • LMDZ5/trunk/libf/phylmd/1D_interp_cases.h

    r2019 r2126  
    7878
    7979      endif ! forcing_toga
     80!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     81! Interpolation DICE forcing
     82!---------------------------------------------------------------------
     83      if (forcing_dice) then
     84
     85       if (prt_level.ge.1) then
     86        print*,'#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_dice=',&
     87     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_dice
     88       endif
     89
     90! time interpolation:
     91      CALL interp_dice_time(daytime,day1,annee_ref                    &
     92     &             ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice     &
     93     &             ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice   &
     94     &             ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice     &
     95     &             ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice &
     96     &             ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof     &
     97     &             ,ustar_prof,psurf_prof,ug_profd,vg_profd           &
     98     &             ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd       &
     99     &             ,omega_profd)
     100!     do l = 1, llm
     101!     print *,'llm l omega_profd',llm,l,omega_profd(l)
     102!     enddo
     103
     104        if (type_ts_forcing.eq.1) ts_cur = tg_prof ! SST used in read_tsurf1d
     105
     106! vertical interpolation:
     107      CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice        &
     108     &         ,th_dice,qv_dice,u_dice,v_dice,o3_dice                   &
     109     &         ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd,omega_profd &
     110     &         ,th_mod,qv_mod,u_mod,v_mod,o3_mod                        &
     111     &         ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
     112!     do l = 1, llm
     113!      print *,'llm l omega_mod',llm,l,omega_mod(l)
     114!     enddo
     115
     116! Les forcages DICE sont donnes /jour et non /seconde !
     117      ht_mod(:)=ht_mod(:)/86400.
     118      hq_mod(:)=hq_mod(:)/86400.
     119      hu_mod(:)=hu_mod(:)/86400.
     120      hv_mod(:)=hv_mod(:)/86400.
     121
     122!calcul de l'advection verticale a partir du omega (repris cas TWPICE, MPL 05082013)
     123!Calcul des gradients verticaux
     124!initialisation
     125      d_t_z(:)=0.
     126      d_q_z(:)=0.
     127      d_u_z(:)=0.
     128      d_v_z(:)=0.
     129      DO l=2,llm-1
     130       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
     131       d_q_z(l)=(q(l+1,1)-q(l-1,1)) /(play(l+1)-play(l-1))
     132       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
     133       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
     134      ENDDO
     135      d_t_z(1)=d_t_z(2)
     136      d_q_z(1)=d_q_z(2)
     137!     d_u_z(1)=u(2)/(play(2)-psurf)/5.
     138!     d_v_z(1)=v(2)/(play(2)-psurf)/5.
     139      d_u_z(1)=0.
     140      d_v_z(1)=0.
     141      d_t_z(llm)=d_t_z(llm-1)
     142      d_q_z(llm)=d_q_z(llm-1)
     143      d_u_z(llm)=d_u_z(llm-1)
     144      d_v_z(llm)=d_v_z(llm-1)
     145
     146!Calcul de l advection verticale:
     147! utiliser omega (Pa/s) et non w (m/s) !! MP 20131108
     148      d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
     149      d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
     150      d_u_dyn_z(:)=omega_mod(:)*d_u_z(:)
     151      d_v_dyn_z(:)=omega_mod(:)*d_v_z(:)
     152
     153! large-scale forcing :
     154!     tsurf = tg_prof    MPL 20130925 commente
     155      psurf = psurf_prof
     156! For this case, fluxes are imposed
     157      fsens=-1*shf_prof
     158      flat=-1*lhf_prof
     159      ust=ustar_prof
     160      tg=tg_prof
     161      print *,'ust= ',ust
     162      do l = 1, llm
     163       ug(l)= ug_profd
     164       vg(l)= vg_profd
     165!       omega(l) = w_prof(l)
     166!      rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
     167!       omega2(l)=-rho(l)*omega(l)
     168!      omega(l) = w_mod(l)*(-rg*rho(l))
     169       omega(l) = omega_mod(l)
     170       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
     171
     172       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
     173       d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
     174       d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
     175       d_u_adv(l) = hu_mod(l)-d_u_dyn_z(l)
     176       d_v_adv(l) = hv_mod(l)-d_v_dyn_z(l)
     177       dt_cooling(l) = 0.0
     178      enddo
     179
     180      endif ! forcing_dice
     181!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     182!---------------------------------------------------------------------
    80183!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    81184!---------------------------------------------------------------------

  • LMDZ5/trunk/libf/phylmd/1D_read_forc_cases.h

    r2117 r2126  
    357357!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    358358!---------------------------------------------------------------------
     359! Forcing from DICE experiment (see file DICE_protocol_vn2-3.pdf)
     360!---------------------------------------------------------------------
     361
     362      if (forcing_dice) then
     363!read DICE forcings
     364      fich_dice='dice_driver.nc'
     365      call read_dice(fich_dice,nlev_dice,nt_dice                    &
     366     &     ,zz_dice,plev_dice,th_dice,qv_dice,u_dice,v_dice,o3_dice &
     367     &     ,shf_dice,lhf_dice,lwup_dice,swup_dice,tg_dice,ustar_dice&
     368     &     ,psurf_dice,ug_dice,vg_dice,ht_dice,hq_dice              &
     369     &     ,hu_dice,hv_dice,w_dice,omega_dice)
     370
     371      write(*,*) 'Forcing DICE lu'
     372
     373!champs initiaux:
     374      do k=1,nlev_dice
     375         th_dicei(k)=th_dice(k)
     376         qv_dicei(k)=qv_dice(k)
     377         u_dicei(k)=u_dice(k)
     378         v_dicei(k)=v_dice(k)
     379         o3_dicei(k)=o3_dice(k)
     380         ht_dicei(k)=ht_dice(k,1)
     381         hq_dicei(k)=hq_dice(k,1)
     382         hu_dicei(k)=hu_dice(k,1)
     383         hv_dicei(k)=hv_dice(k,1)
     384         w_dicei(k)=w_dice(k,1)
     385         omega_dicei(k)=omega_dice(k,1)
     386      enddo   
     387      omega(:)=0.     
     388      omega2(:)=0.
     389      rho(:)=0.
     390! vertical interpolation using TOGA interpolation routine:
     391!      write(*,*)'avant interp vert', t_proftwp
     392!
     393!     CALL interp_dice_time(daytime,day1,annee_ref
     394!    i             ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice
     395!    i             ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice
     396!    i             ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice
     397!    i             ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice
     398!    o             ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof
     399!    o             ,ustar_prof,psurf_prof,ug_profd,vg_profd
     400!    o             ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd
     401!    o             ,omega_profd)
     402
     403      CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice       &
     404     &         ,th_dicei,qv_dicei,u_dicei,v_dicei,o3_dicei             &
     405     &         ,ht_dicei,hq_dicei,hu_dicei,hv_dicei,w_dicei,omega_dicei&
     406     &         ,th_mod,qv_mod,u_mod,v_mod,o3_mod                       &
     407     &         ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
     408
     409! Pour tester les advections horizontales de T et Q, on met w_mod et omega_mod à zero (MPL 20131108)
     410!     w_mod(:,:)=0.
     411!     omega_mod(:,:)=0.
     412
     413!       write(*,*) 'Profil initial forcing DICE interpole',t_mod
     414! Les forcages DICE sont donnes /jour et non /seconde !
     415      ht_mod(:)=ht_mod(:)/86400.
     416      hq_mod(:)=hq_mod(:)/86400.
     417      hu_mod(:)=hu_mod(:)/86400.
     418      hv_mod(:)=hv_mod(:)/86400.
     419
     420! initial and boundary conditions :
     421      write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc
     422      do l = 1, llm
     423! Ligne du dessous à decommenter si on lit theta au lieu de temp
     424       temp(l) = th_mod(l)*(play(l)/pzero)**rkappa
     425!      temp(l) = t_mod(l)
     426       q(l,1) = qv_mod(l)
     427       q(l,2) = 0.0
     428!      print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
     429       u(l) = u_mod(l)
     430       v(l) = v_mod(l)
     431       ug(l)=ug_dice(1)
     432       vg(l)=vg_dice(1)
     433       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
     434!      omega(l) = w_mod(l)*(-rg*rho(l))
     435       omega(l) = omega_mod(l)
     436       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
     437
     438       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
     439!on applique le forcage total au premier pas de temps
     440!attention: signe different de toga
     441       d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
     442!forcage en th
     443!       d_th_adv(l) = ht_mod(l)
     444       d_q_adv(l,1) = hq_mod(l)
     445       d_q_adv(l,2) = 0.0
     446       dt_cooling(l)=0.
     447      enddo     
     448       write(*,*) 'Profil initial forcing DICE interpole temp39',temp(39)
     449     
     450
     451!     ok_flux_surf=.false.
     452      fsens=-1.*shf_dice(1)
     453      flat=-1.*lhf_dice(1)
     454! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par
     455! le coefficient de trainee en surface cd**2=ustar*vent(k=1)
     456! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface
     457! MPL 05082013
     458      ust=ustar_dice(1)
     459      tg=tg_dice(1)
     460      print *,'ust= ',ust
     461      IF (tsurf .LE. 0.) THEN
     462       tsurf= tg_dice(1)
     463      ENDIF
     464      psurf= psurf_dice(1)
     465      solsw_in = (1.-albedo)/albedo*swup_dice(1)
     466      sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1)
     467      PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in
     468      endif !forcing_dice
     469
     470!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     471!---------------------------------------------------------------------
    359472! Forcing from Arm_Cu case                   
    360473! For this case, ifa_armcu.txt contains sensible, latent heat fluxes

  • LMDZ5/trunk/libf/phylmd/clcdrag.F90

    r1907 r2126  
    9292        IF (.NOT.zxli) THEN
    9393           zscf = SQRT(1.+cd*ABS(zri(i)))
    94            FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1)
     94           FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min)
    9595           zcfm1(i) = zcdn(i) * FRIV
    96            FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 )
     96           FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), f_ri_cd_min )
    9797!!$  PB          zcfh1(i) = zcdn(i) * FRIH
    9898!!$ PB           zcfh1(i) = f_cdrag_stable * zcdn(i) * FRIH
    … …  
    110110           zucf = 1./(1.+3.0*cb*cc*zcdn(i)*SQRT(ABS(zri(i)) &
    111111                *(1.0+zgeop1(i)/(RG*rugos(i)))))
    112            zcfm2(i) = zcdn(i)*amax1((1.-2.0*cb*zri(i)*zucf),0.1)
    113 !!$PB            zcfh2(i) = zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
    114            zcfh2(i) = f_cdrag_ter*zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
     112           zcfm2(i) = zcdn(i)*amax1((1.-2.0*cb*zri(i)*zucf),f_ri_cd_min)
     113!!$PB            zcfh2(i) = zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
     114           zcfh2(i) = f_cdrag_ter*zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
    115115           pcfm(i) = zcfm2(i)
    116116           pcfh(i) = zcfh2(i)

  • LMDZ5/trunk/libf/phylmd/clesphys.h

    r2114 r2126  
    3434       REAL cdmmax, cdhmax
    3535!IM param. stabilite s/ terres et en dehors
    36        REAL ksta, ksta_ter
     36       REAL ksta, ksta_ter, f_ri_cd_min
    3737!IM ok_kzmin : clef calcul Kzmin dans la CL de surface cf FH
    3838       LOGICAL ok_kzmin
    … …  
    8787     &     , CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt                     &
    8888     &     , CH4_ppb_per, N2O_ppb_per, CFC11_ppt_per, CFC12_ppt_per     &
    89      &     , cdmmax, cdhmax, ksta, ksta_ter                             &
     89     &     , cdmmax, cdhmax, ksta, ksta_ter, f_ri_cd_min                &
    9090     &     , fmagic, pmagic                                             &
    9191     &     , f_cdrag_ter,f_cdrag_oce,f_rugoro                           &

  • LMDZ5/trunk/libf/phylmd/coefcdrag.F90

    r2011 r2126  
    6363      REAL, dimension(klon) :: zcfm2, zcfh2
    6464      REAL, dimension(klon) :: trm0, trm1
     65
    6566!-------------------------------------------------------------------------
    6667      REAL :: fsta, fins, x
    … …  
    110111         IF (.NOT.zxli) THEN
    111112           zscf(i) = SQRT(1.+CD*ABS(zri1(i)))
    112            friv(i) = max(1. / (1.+2.*CB*zri1(i)/ zscf(i)), 0.1)
     113           friv(i) = max(1. / (1.+2.*CB*zri1(i)/ zscf(i)), f_ri_cd_min)
    113114           zcfm1(i) = cdran(i) * friv(i)
    114            frih(i) = max(1./ (1.+3.*CB*zri1(i)*zscf(i)), 0.1 )
     115           frih(i) = max(1./ (1.+3.*CB*zri1(i)*zscf(i)), f_ri_cd_min )
    115116!           zcfh1(i) = cdran(i) * frih(i)
    116117           zcfh1(i) = f_cdrag_ter*cdran(i) * frih(i)
    … …  
    130131           zucf(i) = 1./(1.+3.0*CB*CC*cdran(i)*SQRT(ABS(zri1(i)) &
    131132                 *(1.0+zdphi(i)/(RG*rugos(i)))))
    132            zcfm2(i) = cdran(i)*max((1.-2.0*CB*zri1(i)*zucf(i)),0.1)
    133 !           zcfh2(i) = cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),0.1)
    134            zcfh2(i) = f_cdrag_ter*cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),0.1)
     133           zcfm2(i) = cdran(i)*max((1.-2.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
     134!           zcfh2(i) = cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
     135           zcfh2(i) = f_cdrag_ter*cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),f_ri_cd_min)
    135136           cdram(i) = zcfm2(i)
    136137           cdrah(i) = zcfh2(i)

  • LMDZ5/trunk/libf/phylmd/conf_phys_m.F90

    r2114 r2126  
    160160    REAL,SAVE :: solarlong0_omp
    161161    INTEGER,SAVE :: top_height_omp,overlap_omp
    162     REAL,SAVE :: cdmmax_omp,cdhmax_omp,ksta_omp,ksta_ter_omp
     162    REAL,SAVE :: cdmmax_omp,cdhmax_omp,ksta_omp,ksta_ter_omp,f_ri_cd_min_omp
    163163    LOGICAL,SAVE :: ok_kzmin_omp
    164164    REAL, SAVE ::  fmagic_omp, pmagic_omp
    … …  
    11421142    ksta_ter_omp = 1.0e-10
    11431143    call getin('ksta_ter',ksta_ter_omp)
     1144
     1145    !Config Key  = f_ri_cd_min
     1146    !Config Desc =
     1147    !Config Def  = 0.1
     1148    !Config Help =
     1149    !
     1150    f_ri_cd_min_omp = 0.1
     1151    call getin('f_ri_cd_min',f_ri_cd_min_omp)
    11441152
    11451153    !
    … …  
    18311839    ksta = ksta_omp
    18321840    ksta_ter = ksta_ter_omp
     1841    f_ri_cd_min = f_ri_cd_min_omp
    18331842    ok_kzmin = ok_kzmin_omp
    18341843    fmagic = fmagic_omp
    … …  
    20662075    write(lunout,*)' ksta = ',ksta
    20672076    write(lunout,*)' ksta_ter = ',ksta_ter
     2077    write(lunout,*)' f_ri_cd_min = ',f_ri_cd_min
    20682078    write(lunout,*)' ok_kzmin = ',ok_kzmin
    20692079    write(lunout,*)' fmagic = ',fmagic

  • LMDZ5/trunk/libf/phylmd/flux_arp.h

    r1907 r2126  
    33!
    44      logical :: ok_flux_surf
     5      logical :: ok_prescr_ust !for prescribed ustar
    56      real :: fsens
    67      real :: flat
     8      real :: ust
     9      real :: tg
    710
    8       common /flux_arp/fsens,flat,ok_flux_surf
     11      common /flux_arp/fsens,flat,ust,tg,ok_flux_surf,ok_prescr_ust
    912
    1013!$OMP THREADPRIVATE(/flux_arp/)

  • LMDZ5/trunk/libf/phylmd/lmdz1d.F90

    r2117 r2126  
    112112        logical :: forcing_twpice  = .false.
    113113        logical :: forcing_amma    = .false.
     114        logical :: forcing_dice    = .false.
    114115        logical :: forcing_GCM2SCM = .false.
    115116        logical :: forcing_GCSSold = .false.
    … …  
    178179!---------------------------------------------------------------------
    179180      logical :: ok_writedem =.true.
     181      real :: sollw_in = 0.
     182      real :: solsw_in = 0.
    180183     
    181184!---------------------------------------------------------------------
    … …  
    262265!             initial profiles from AMMA nc file
    263266!             LS convergence, omega and surface fluxes imposed from AMMA file 
     267!forcing_type = 7 ==> forcing_dice = .true.
     268!             initial profiles and large scale forcings in dice_driver.nc
     269!             Different stages: soil model alone, atm. model alone
     270!             then both models coupled
    264271!forcing_type = 40 ==> forcing_GCSSold = .true.
    265272!             initial profile from GCSS file
    … …  
    296303      elseif (forcing_type .eq.6) THEN
    297304       forcing_amma = .true.
     305      elseif (forcing_type .eq.7) THEN
     306       forcing_dice = .true.
    298307      elseif (forcing_type .eq.40) THEN
    299308       forcing_GCSSold = .true.
    … …  
    319328
    320329        type_ts_forcing = 0
    321         if (forcing_toga.or.forcing_sandu.or.forcing_astex)                 &
     330        if (forcing_toga.or.forcing_sandu.or.forcing_astex .or. forcing_dice)                 &
    322331     &    type_ts_forcing = 1
    323332
    … …  
    366375! Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
    367376      IF(forcing_type .EQ. 6) fnday=64800./86400.
     377!     IF(forcing_type .EQ. 6) fnday=50400./86400.
    368378      annee_ref = anneeref
    369379      mois = 1
    … …  
    391401     & (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma              &
    392402     & ,day_ju_ini_amma)
    393 
     403      ELSEIF (forcing_type .eq.7) THEN
     404! Convert the initial date of DICE to Julian day
     405      call ymds2ju                                                         &
     406     & (year_ini_dice,mth_ini_dice,day_ini_dice,heure_ini_dice             &
     407     & ,day_ju_ini_dice)
    394408      ELSEIF (forcing_type .eq.59) THEN
    395409! Convert the initial date of Sandu case to Julian day
    … …  
    405419
    406420      ELSEIF (forcing_type .eq.61) THEN
    407 
    408421! Convert the initial date of Arm_cu case to Julian day
    409422      call ymds2ju                                                          &

  • LMDZ5/trunk/libf/phylmd/pbl_surface_mod.F90

    r1932 r2126  
    319319    REAL, DIMENSION(klon),        INTENT(OUT)       :: alb2_m     ! mean albedo in near IR SW interval
    320320    ! Martin
    321         REAL, DIMENSION(klon),        INTENT(OUT)       :: alb3_lic
     321    REAL, DIMENSION(klon),        INTENT(OUT)       :: alb3_lic
    322322    ! Martin
    323323    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens     ! sensible heat flux at surface with inversed sign
    … …  
    475475    LOGICAL, PARAMETER                 :: check=.FALSE.
    476476    REAL, DIMENSION(klon)              :: Kech_h       ! Coefficient d'echange pour l'energie
     477    REAL                               :: vent
    477478
    478479! For debugging with IOIPSL
    … …  
    797798            yts, yqsurf, yrugos, &
    798799            ycdragm, ycdragh )
     800! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
     801     IF (ok_prescr_ust) then
     802      DO i = 1, knon
     803       print *,'ycdragm avant=',ycdragm(i)
     804       vent= sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))
     805!      ycdragm(i) = ust*ust/(1.+(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1)))
     806!      ycdragm(i) = ust*ust/((1.+sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))) &
     807!     *sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1)))
     808       ycdragm(i) = ust*ust/(1.+vent)/vent
     809       print *,'ycdragm ust yu yv apres=',ycdragm(i),ust,yu(i,1),yv(i,1)
     810      ENDDO
     811     ENDIF
     812
    799813
    800814!****************************************************************************************
    … …  
    905919               y_flux_u1, y_flux_v1 )
    906920               
     921! Special DICE MPL 05082013
     922       IF (ok_prescr_ust) THEN
     923!         ysnow(:)=0.
     924!         yqsol(:)=0.
     925!         yagesno(:)=50.
     926!         ytsoil(:,:)=300.
     927!         yz0_new(:)=0.001
     928!         yalb1_new(:)=0.22
     929!         yalb2_new(:)=0.22
     930!         yevap(:)=flat/RLVTT
     931!         yfluxlat(:)=-flat
     932!         yfluxsens(:)=-fsens
     933!         yqsurf(:)=0.
     934!         ytsurf_new(:)=tg
     935!         y_dflux_t(:)=0.
     936!         y_dflux_q(:)=0.
     937          y_flux_u1(:)=ycdragm(:)*(1.+sqrt(yu(:,1)*yu(:,1)+yv(:,1)*yv(:,1)))*yu(:,1)*ypplay(:,1)/RD/yt(:,1)
     938          y_flux_v1(:)=ycdragm(:)*(1.+sqrt(yu(:,1)*yu(:,1)+yv(:,1)*yv(:,1)))*yv(:,1)*ypplay(:,1)/RD/yt(:,1)
     939      ENDIF
     940
    907941     
    908942       CASE(is_lic)

  • LMDZ5/trunk/libf/phylmd/stdlevvar.F90

    r1907 r2126  
    5353      REAL, dimension(klon), intent(out) :: u_10m, t_10m, q_10m
    5454!-------------------------------------------------------------------------
     55      include "flux_arp.h"
    5556      include "YOMCST.h"
    5657!IM PLUS
    … …  
    101102 &                   ts1, qsurf, rugos, okri, ri1,  &         
    102103 &                   cdram, cdrah, cdran, zri1, pref)           
     104! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
     105     IF (ok_prescr_ust) then
     106      DO i = 1, knon
     107       print *,'cdram avant=',cdram(i)
     108       cdram(i) = ust*ust/speed(i)/speed(i)
     109       print *,'cdram ust speed apres=',cdram(i),ust,speed
     110      ENDDO
     111     ENDIF
    103112!
    104113!---------Star variables----------------------------------------------------

  • LMDZ5/trunk/libf/phylmd/surf_land_orchidee_mod.F90

    r2011 r2126  
    614614          correspond(ij,jj) = igrid
    615615       ENDDO
     616!sonia : Les mailles des voisines doivent etre toutes egales (pour couplage orchidee)
     617       IF (knon_glo == 1) THEN
     618         igrid = 1
     619         DO i = 1,8
     620           neighbours_glo(igrid, i) = igrid
     621         ENDDO
     622       ELSE
     623       print*,'sonia : knon_glo,ij,jj', knon_glo, ij,jj
    616624       
    617625       DO igrid = 1, knon_glo
    … …  
    636644          ENDDO
    637645       ENDDO
     646       ENDIF !fin knon_glo == 1
    638647
    639648    ENDIF
Note: See TracChangeset for help on using the changeset viewer.