source: LMDZ6/branches/LMDZ_DECOUPLE/libf/phymar/PHY_SISVAT_INI.f90 @ 5434

Last change on this file since 5434 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR


Integration of MAR physics

File size: 18.8 KB
Line 
1      subroutine PHY_SISVAT_INI
2
3!------------------------------------------------------------------------------+
4!                                                         Sat 29-Jun-2013  MAR |
5!     SubRoutine PHY_SISVAT_INI initializes                                    |
6!                    SISVAT    (Soil                                           |
7!                               Ice                                            |
8!                               Snow                                           |
9!                               Vegetation                                     |
10!                               Atmosphere                                     |
11!                               Transfer   Scheme)                             |
12!                                                                              |
13!     version 3.p.4.1 created by H. Gallee,               Mon  4-Feb-2013      |
14!           Last Modification by H. Gallee,               Sat 29-Jun-2013      |
15!                                                                              |
16!------------------------------------------------------------------------------+
17
18
19! Global  Variables
20! =================
21
22      use Mod_Real
23      use Mod_PHY____grd
24      use Mod_PHY____dat
25      use Mod_PHY____kkl
26      use Mod_PHY_DY_kkl
27      use Mod_SISVAT_ctr
28      use Mod_SISVAT_grd
29      use Mod_SISVAT_dat
30      use Mod_SISVAT_loc
31      use Mod_SISVAT_kkl
32      use Mod_SISVAT_gpt
33
34
35
36      IMPLICIT NONE
37
38
39      integer                 ::  isl   ,n     ,k     ,l     ,jo    ,no
40      integer                 ::  i     ,j     ,ikp   ,nm
41
42
43
44! Normalized Decay of the Surficial Water Content: data
45! ===============================================
46
47      real(kind=real8)  ::  c1_zuo = 12.960e+4      !  Run Off Parameters
48      real(kind=real8)  ::  c2_zuo =  2.160e+6      !  86400*1.5 day     ...*25 days
49                                                    !  86400*0.3 day    (Modif. ETH Camp)
50      real(kind=real8)  ::  c3_zuo =  1.400e+2      !  (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)
51
52!     real(kind=real8)  ::  c1_zuo =  2.796e+4      !  Run Off Parameters (Tuning)
53
54
55
56 
57! INITIALISATION:  BEGIN  ++++++++++++++++++++++++++++++++++++++++++++ 
58
59
60! Initialisation of Mod_SISVAT_dat (Ocean Surface Parameters)
61! ================================
62
63! #SI VarSST =   0.
64! #OP VarSST =   1.                  ! Variable (0.) / Fixed    (1.) SSTs
65      FixSST =   1.-VarSST           ! Fixed    (1.) / Variable (0.) SSTs
66      SSTnud = exp(-dt__SV/2.592e6)  ! SST    Nudging:
67!                                    ! e-folding time: 30 Days
68
69! #SI TOF_SV = TocnSI                ! Ocn Grid Cell Freez.Temperature  [K]
70! #RE TOF_SV = 271.35 + eps6         !
71
72
73 
74
75! Initialization of Mod_SISVAT_grd
76! ================================
77
78! Parameters used in the Interpolation of V(10 m)
79! -----------------------------------------------
80
81        IF   (hsigma(1).GT.10.          )                           THEN
82          k = 0
83  301     CONTINUE
84          k = k + 1
85            !if (hsigma(k).LT.10.OR.k.gt.mzp)                   GO TO 300
86            if (hsigma(k).LT.10.OR.k.ge.mzp)                   GO TO 300
87                                                               GO TO 301
88  300     CONTINUE
89          k_SL = k
90
91          IF (k_SL.EQ.mzp)                                          THEN
92              r_SL10 = log(10.           / 0.002)   &! 0.002: typical Z0
93     &                /log(hsigma(k_SL)  / 0.002)    !
94
95          ELSE
96              r_SL10 =    (10.           - hsigma(k_SL))               &
97     &                /   (hsigma(k_SL-1)- hsigma(k_SL))
98          END IF
99        ELSE
100              k_SL   =     mzp
101              r_SL10 =      1.
102        END IF
103
104
105
106! Level of negligible blown Particles Concentration ( z_zb ~ 25 magl)      SNOW
107! -------------------------------------------------                   .OR. DUST
108
109! #AE        k_zb  =mzp
110! #AE 11 CONTINUE
111! #AE    IF (hsigma(k_zb  ).GT.z_zb.OR.k_zb  .LE.1)             GO TO 10
112! #AE        k_zb  =k_zb  -1
113! #AE                                                           GO TO 11
114! #AE 10 CONTINUE
115! #AE        write(6,1000)             k_zb
116! #AE 1000   format(/,' BS : Level of negligible '                     &
117! #AE&                     ,'blown Particles Concentration is',i4      &
118! #AE&                     ,' (i.e., ~  25. magl)',/)
119
120
121
122
123! Initialization of SISVAT constants and parameters
124! =================================================
125
126!           **********
127      CALL  SISVAT_ini
128!           **********
129
130
131
132
133! Initialization of Mod_SISVAT_kkl
134! ================================
135
136! Surface Fall Line Slope
137! -----------------------
138
139        IF (kcolp .EQ. 1)                                           THEN
140            slopAP    (1)    =      slop1d
141          DO n   = 1,mwp
142            slopSV(1  ,n)    =      slop1d
143          ENDDO
144        ELSE
145          DO ikp=1,kcolp
146          DO nm =1,mwp
147            slopSV(ikp,nm)   = slopAP(ikp)
148          END DO
149          END DO
150        END IF
151
152        IF (SnoMod)                                                 THEN
153          DO ikp=1,kcolp
154          DO nm =1,mwp
155            SWf_SV(ikp,nm)      =                            &! Normalized Decay of the
156     &       exp(-dt__SV                                  &! Surficial Water Content
157     &          /(c1_zuo                                  &! Zuo and Oerlemans 1996,
158     &           +c2_zuo*exp(-c3_zuo*slopSV(ikp,nm) )))       ! J.Glacio. 42, 305--317
159          END DO
160          END DO
161        END IF
162
163
164
165
166
167! OUTPUT point (i,j,n) coordinates
168! --------------------------------
169
170! stdout
171! ~~~~~~
172! Martin CONTROL
173          iwr_SV = 23
174          jwr_SV = 17
175          nwr_SV = 1
176          !iwr_SV = 1
177          !jwr_SV = 1
178          !nwr_SV = 1
179! Martin CONTROL
180
181! Particular txt files (for a regular spacing in 1D arrays)
182! ~~~~~~~~~~~~~~~~~~~~
183        DO  ikp=1,kcolp
184        DO  nm =1,mwp
185          lwriSV(ikp,nm) = 0
186        END DO
187        END DO
188
189          jo             = kcolp*mwp/nbwri
190          jo             = max(1,jo)
191          no             =               0
192          ikp            =      -jo /    2
193 11     CONTINUE
194          nm             =       1
195          no             = no  + 1
196          ikp            = ikp + jo
197          ikp            = max(    1,ikp)
198          ikp            = min(kcolp,ikp)
199
200 101    CONTINUE
201          i___SV(no)     = ii__AP(ikp)
202          j___SV(no)     = jj__AP(ikp)
203          n___SV(no)     = nm
204          lwriSV(ikp,nm) = 1
205          nm             = nm  + 1
206          no             = no  + 1
207        IF (nm.LE.mwp .AND. no.LE.nbwri)                       GO TO 101
208          no          = no  - 1
209
210        IF (no.LT.nbwri)                                       GO TO 11
211
212
213
214
215! +++++  ++++++++++++++  +++++  ++++++++++++++++++++++++++++++++++++++++
216! FIRST  INITIALISATION: BEGIN
217! +++++  ++++++++++++++  +++++
218
219!                 ===========
220        IF       (it_EXP.EQ.1)                                      THEN
221!                 ===========
222
223
224! Ocean                          1st Initialization (Grid Cells)
225! -------------------------------------------------
226
227! #SI     DO ikp=1,kcolp
228! #SI       IF  (MaskSV_gpt(ikp) .EQ. 0 .AND. FracSV(ikp,1) .lt. 1.)  THEN
229! #SI            i   = ii__AP(ikp)
230! #SI            j   = jj__AP(ikp)
231! #SI          DO n=1,mwp
232! #SI            write(6,6000)i,j,n,FracSV(ikp,nm)
233! #SI 6000       format(' WARNING on Grid Point',2i4,' Mosaic',i3      &
234! #SI&                 ,' Fraction = f4.3,' : ISLAND excluded')
235! #SI            FracSV    (ikp,nm)  =   0.
236! #SI            iVgTSV    (ikp,nm)  =   0
237! #SI          END DO
238! #SI            FracSV    (ikp,1)   =   1.
239! #SI       END IF
240! #SI     END DO
241
242
243! Prescription from SST
244! ---------------------
245
246! #SI     DO ikp=1,kcolp
247! #SI        FraOcn          =(TsisSV(ikp,1,0)-TOF_SV)/TSIdSV           ! Open Ocean
248! #IP        FraOcn          = 1.0000         -sif_SB(ikp)              ! Prescribed
249! #SI        FraOcn          = min(  un_1,FraOcn)                       !      Fract.
250! #SI        FraOcn          = max(OcnMin,FraOcn)                       !
251! #SI        FracSV(ikp,1)   = LSMask(ikp,1)     *    FracSV(ikp,1)    &! New  Ocean
252! #SI&                     +(1-LSMask(ikp,1)    )*    FraOcn            !
253! #SI        SrfSIC          =                        FracSV(ikp,2)     ! Old  Sea Ice
254! #SI        SIc0OK          = max(zer0, sign(un_1,   SrfSIC-eps6))     !
255! #SI        FracSV(ikp,2)   = LSMask(ikp,1)     *    FracSV(ikp,2)    &! New  Sea Ice
256! #SI&                     +(1-LSMask(ikp,1)    )*(1.-FraOcn)           !
257! #SI        SIceOK          = max(zer0, sign(un_1,   FracSV(ikp,2)    &!
258! #SI&                                                      -eps6))     !
259
260! Sea-Ice Vertical Discretization
261! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
262! #SI        isnoSV    (ikp,2)     =                                   &
263! #SI&       isnoSV    (ikp,2)                   *   LSMask(ikp,2)     &
264! #SI&     +(isnoSV    (ikp,2)  * SIc0OK                               &
265! #SI&      +     3          *(1.-SIc0OK)*SIceOK)*(1-LSMask(ikp,2) )
266
267! #SI        iiceSV    (ikp,2)     =                                   &
268! #SI&       iiceSV    (ikp,2)                   *   LSMask(ikp,2)     &
269! #SI&     +(iiceSV    (ikp,2)     *SIc0OK                             &
270! #SI&      +     3          *(1.-SIc0OK)*SIceOK)*(1-LSMask(ikp,2) )
271! #SI        ispiSV    (ikp,2)     =                 iiceSV(ikp,2) 
272
273! #SI       DO l=1,nsnow
274! #SI        dzsnSV    (ikp,2,l) =                                     &
275! #SI&       dzsnSV    (ikp,2,l)                 *   LSMask(ikp,2)     &
276! #SI&     +(dzsnSV    (ikp,2,l) *SIc0OK                               &
277! #SI&      +dzSIce(min(4,l))*(1.-SIc0OK)*SIceOK)*(1-LSMask(ikp,2) )
278
279! #SI        TsisSV    (ikp,2,l) =                                     &
280! #SI&       TsisSV    (ikp,2,l)                 *   LSMask(ikp,2)     &
281! #SI&     +(TsisSV    (ikp,2,l) *    SIc0OK                           &
282! #SI&      +TsisSV    (ikp,1,0) *(1.-SIc0OK)   )*(1-LSMask(ikp,2) )
283
284! #SI        ro__SV    (ikp,2,l) =                                     &
285! #SI&       ro__SV    (ikp,2,l)                 *   LSMask(ikp,2)     &
286! #SI&     +(ro__SV    (ikp,2,l) *SIc0OK                               &
287! #SI&      +rhoIce          *(1.-SIc0OK)*SIceOK)*(1-LSMask(ikp,2) )
288
289! #SI        G1snSV    (ikp,2,l) =                                     &
290! #SI&       G1snSV    (ikp,2,l)                 *   LSMask(ikp,2)     &
291! #SI&     +(G1snSV    (ikp,2,l) *SIc0OK                               &
292! #SI&      +G1_dSV          *(1.-SIc0OK)*SIceOK)*(1-LSMask(ikp,2) )
293
294! #SI        G2snSV    (ikp,2,l) =                                     &
295! #SI&       G2snSV    (ikp,2,l)                 *   LSMask(ikp,2)     &
296! #SI&     +(G2snSV    (ikp,2,l) *SIc0OK                               &
297! #SI&      +30.             *(1.-SIc0OK)*SIceOK)*(1-LSMask(ikp,2) )
298
299! #SI        istoSV    (ikp,2,l) =                                     &
300! #SI&       istoSV    (ikp,2,l)                 *   LSMask(ikp,2)     &
301!!#SI&new  +(istoSV    (ikp,2,l) *SIc0OK                               &
302!!#SI&new   +istdSV(2)       *(1.-SIc0OK)*SIceOK)*(1-LSMask(ikp,2) )   &
303! #SI&     + istdSV(2)           *                (1-LSMask(ikp,2) )
304! #SI       END DO
305! #SI       DO l=-nsoil,0
306! #SI        TsisSV    (ikp,2,l) =                                     &
307! #SI&       TsisSV    (ikp,2,l)                 *   LSMask(ikp,2)     &
308! #SI&     +(TsisSV    (ikp,2,l) *    SIc0OK                           &
309! #SI&      +TsisSV    (ikp,1,l) *(1.-SIc0OK)   )*(1-LSMask(ikp,2) )
310
311! #SI        eta_SV    (ikp,2,l) =                                     &
312! #SI&       eta_SV    (ikp,2,l)                 *   LSMask(ikp,2)     &
313! #SI&     + eta_SV    (ikp,2,l) *    SIc0OK     *(1-LSMask(ikp,2) )
314!                                 No Pore in Ice => No Water
315! #SI       END DO
316
317! OUTPUT of Sea Ice Fraction, Temperature, Discretization
318! #si        write(6,6001) Day_TU,LabMon(Mon_TU),YearTU                &
319! #si&                    ,HourTU,minuTU,Sec_TU ,TsisSV(ikp,1,0)       &
320! #si&                    ,FraOcn,FracSV(ikp,1) ,TsisSV(ikp,2,0)       &
321! #si&                    ,iiceSV(ikp,2)        ,isnoSV(ikp,2) 
322! #si 6001        format(/,98('_')                                     &
323! #si&         ,/, i3,'-',a3,'-',i4,3(':',i2)                          &
324! #si&         ,   2x,'T OCN = ',f7.3,4x,'% OCN = ',f7.3,'(',f4.3,')'  &
325! #si&         ,   2x,'T ICE = ',f7.3                                  &
326! #si&         ,/,43x,'NbIce = ',i3, 11x,'NbSno = ',i3)     
327
328! #SI     END DO
329
330
331! SBL                            1st Initialization
332! -------------------------------------------------
333
334! Influence of the Angle(Wind,Sastrugi) (Andreas, 1995, CCREL report 95-16)
335! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
336! #Za     DO ikp=1,kcolp
337! #Za     DO nm =1,mwp
338! #Za       ua0_SV(ikp,nm)  = ua__SV(ikp,nm)
339! #Za       va0_SV(ikp,nm)  = va__SV(ikp,nm)
340! #Za     END DO
341! #Za     END DO
342
343! H2O  Upward IR Flux                Initialization
344! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
345          DO ikp=1,kcolp
346            WE2aSV_gpt(ikp) = 0.0
347          END DO
348
349
350! GUESS, eventually from DATA   (1st Initialization)
351! --------------------------------------------------
352
353          DO ikp = 1,kcolp
354          DO nm  = 1,mwp
355
356            IRs_SV(ikp,nm)    =-StefBo*Ta__DY(ikp,mzpp) *Ta__DY(ikp,mzpp)     &! Upward IR Flux
357     &                                *Ta__DY(ikp,mzpp) *Ta__DY(ikp,mzpp)      !
358
359
360! SBL                            1st Initialization (Mosaics)
361! -------------------------------------------------
362
363! Drag Coefficient               1st Initialization
364! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
365            rCDmSV(ikp,nm)    =    0.04
366            rCDhSV(ikp,nm)    =    0.04
367
368! Turbulent Scales               1st Initialization
369! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
370            us__SV(ikp,nm)    = rCDmSV(ikp,nm) * WindSV(ikp,nm,mzp)
371            uts_SV(ikp,nm)    = rCDhSV(ikp,nm) *(Ta__DY(ikp,mzp)-Ta__DY(ikp,mzpp) ) &
372     &                        * us__SV(ikp,nm)
373            uqs_SV(ikp,nm)    =    0.0
374            uqs_SV_gpt(ikp)   =    0.0                                   !  redondance sur n, sans consequences
375            uss_SV(ikp,nm)    =    0.0
376! #BS       ussxSV(ikp,nm)    =    0.0
377
378! Orography Roughness Lengths
379! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
380! #ZO       Z0roSV(ikp,nm)    = min(z0__SV(ikp,nm),hsigma(mzp)/3.)
381
382! Roughness Lengths              1st Initialization
383! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
384! #ZS       Z0SaSV(ikp,nm)    =    0.0                                   !  z0(Sastrugi  h)
385! #ZM      DO k=1,ntavz
386! #ZM       z0_mem(ikp,nm,k)  =    0.0001
387! #ZM       r0_mem(ikp,nm,k)  =    0.0001
388! #ZM       b0_mem(ikp,nm,k)  =    0.0001
389! #ZM      END DO
390
391! SBL Wind Speed and Vertical Temperature Gradient
392! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
393! #AA      DO k=1,ntave
394! #AW       V__mem(ikp,nm,k)  =    WindSV(ikp,nm,mzp)
395! #AH       T__mem(ikp,nm,k)  =    0.0000
396! #AA      END DO
397
398! SBL Water Vapor Flux Limitor
399! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
400! #WL      DO n =1,nLimi
401! #WL        WL_mem(ikp,nm,n)  =   1.
402! #WL      END DO
403
404
405! Frazil                         1st Initialization
406! -------------------------------------------------
407
408             HFraSV(ikp,nm)    =   0.
409
410
411! Roughness Length               1st Initialization
412! -------------------------------------------------
413
414             Z0roSV(ikp,nm)    =   z0__SV(ikp,nm) 
415
416
417! Vegetation                     1st Initialization
418! -------------------------------------------------
419
420             TvegSV(ikp,nm)    = Ta__DY(ikp,mzpp)     ! Vegetation skin Temperature
421             snCaSV(ikp,nm)    =          0.          ! Canopy intercepted Snow
422             rrCaSV(ikp,nm)    =          0.          ! Canopy intercepted Raiw
423             psivSV(ikp,nm)    =          0.          ! Leaf   Water    Potential     [m]
424
425
426! Blowing Snow                   1st Initialization
427! -------------------------------------------------
428
429! #Za        VVs_SV(ikp,nm)    =         10.          ! Wind Speed,         (Sastrugi), Relevance  [m/s]
430! #Za        RRsxSV(ikp,nm)    =          1.          ! Wind Speed Counter, (Sastrugi), Relevance    [-]
431! #Za        DDsxSV(ikp,nm)    =          0.          ! Wind Direction    , (Sastrugi), Relevance    [-]
432
433             BufsSV(ikp,nm)    =          0.          ! Fallen Snow Buffer          (Mosaic)   [mm w.e.]
434             ussbSV(ikp,nm)    =          0.          ! Eroded Snow Buffer          (Mosaic)   [mm w.e.]
435
436
437
438! Snow Pack                      1st Initialization
439! -------------------------------------------------
440
441             zWE_SV(ikp,nm)    =            0.        ! Thickness
442             qSalSV(ikp,nm)    =            0.        ! Saltating Particles Concentration
443             BrosSV(ikp,nm)    =          300.        ! Buffer Snow Layer: initial density: Polar Snow
444             BG1sSV(ikp,nm)    =       G1_dSV         ! Buffer Snow Layer: initial G1     : Polar Snow
445             BG2sSV(ikp,nm)    =       ADSdSV         ! Buffer Snow Layer: initial G2     : Polar Snow
446                                                   ! Buffer Snow Layer  initial characteristics for 0-mass
447
448             wes_SV(ikp,nm)    =            0.        ! Depo. / Subli.
449             wem_SV(ikp,nm)    =            0.        ! Melting
450             wer_SV(ikp,nm)    =            0.        ! Refreezing
451             wee_SV(ikp,nm)    =            0.        ! Evapotranspiration
452
453
454! Surficial  Water                   Initialization
455! -------------------------------------------------
456
457             SWf_SV(ikp,nm)    =            0.        ! Normalized Decay of Surficial Water Content  [-]
458             rusnSV(ikp,nm)    =            0.        ! Surficial Water Mass        (Mosaic)   [mm w.e.]
459             SWS_SV(ikp,nm)    =            1.        ! Surficial Water Status      (Mosaic)   [0,1=m,f]
460                                                   ! Freezng Initial Status   is assumed
461 
462
463! Cumulative Run-Off                 Initialization
464! -------------------------------------------------
465
466             ROF_SV(ikp,nm)       =  0.
467
468
469! Soil Roots                     1st Initialization
470! -------------------------------------------------
471
472           DO isl = -nsoil,0
473             Rootsv(ikp,nm,isl) =    0.
474           END DO
475
476          END DO
477          END DO
478
479!                 ===========
480        END IF ! (it_EXP.EQ.0)
481!                 ===========
482!
483! +++++  ++++++++++++++  +++
484! FIRST  INITIALISATION: END
485! +++++  ++++++++++++++  +++  ++++++++++++++++++++++++++++++++++++++++++
486
487
488
489
490! OUTPUT Files Set-Up
491! =======================
492
493! #v0   open(unit=50,status='unknown',file='PHY_SISVAT.v0')
494! #v0   rewind    50
495
496
497
498
499! OUTPUT
500! ======
501
502      IF (kcolp.EQ.1)                                               THEN
503         write(6,6) (iVgTSV    (1,nm),nm=1,mwp)
504 6       format(/ ,' Vegetation Type: ',6i9)
505      END IF
506
507
508
509! +++ INITIALISATION:  END  ++++++++++++++++++++++++++++++++++++++++++ 
510
511      return
512      end subroutine PHY_SISVAT_INI
Note: See TracBrowser for help on using the repository browser.