source: LMDZ.3.3/trunk/libf/phylmd/interface_surf.F90 @ 191

Last change on this file since 191 was 175, checked in by lmdz, 24 years ago

Rajout de l'include indicesol.inc F90
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.7 KB
Line 
1
2  MODULE interface_surf
3
4! Ce module regroupe toutes les routines gerant l'interface entre le modele
5! atmospherique et les modeles de surface (sols continentaux, oceans, glaces)
6! Les routines sont les suivantes:
7!
8!   interfsurf_*: routines d'aiguillage vers les interfaces avec les
9!                 differents modeles de surface
10!   interfsol\
11!             > routines d'interface proprement dite
12!   interfoce/
13!
14!   interfstart: routine d'initialisation et de lecture de l'etat initial
15!                "interface"
16!   interffin  : routine d'ecriture de l'etat de redemmarage de l'interface
17!
18!
19! L. Fairhead, LMD, 02/2000
20
21  USE ioipsl
22!  USE constantes
23
24  IMPLICIT none
25
26  PRIVATE
27  PUBLIC :: interfsurf,interfsurf_hq
28
29  INTERFACE interfsurf
30    module procedure interfsurf_hq, interfsurf_vent
31  END INTERFACE
32
33  INTERFACE interfoce
34    module procedure interfoce_cpl, interfoce_slab, interfoce_lim
35  END INTERFACE
36
37#include "YOMCST.inc"
38#include "indicesol.inc"
39
40
41! run_off      ruissellement total
42  real, allocatable, dimension(:),save    :: run_off
43
44
45  CONTAINS
46!
47!############################################################################
48!
49  SUBROUTINE interfsurf_hq(itime, dtime, jour, &
50      & klon, nisurf, knon, knindex, rlon, rlat, &
51      & debut, lafin, ok_veget, &
52      & zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, &
53      & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
54      & precip_rain, precip_snow, lwdown, swnet, swdown, &
55      & tsurf, p1lay, cal, beta, coef1lay, ps, radsol, dif_grnd, &
56      & ocean, &
57      & evap, fluxsens, fluxlat, dflux_l, dflux_s, &             
58      & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new)
59
60
61! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general
62! (sols continentaux, oceans, glaces) pour les fluxs de chaleur et d'humidite.
63! En pratique l'interface se fait entre la couche limite du modele
64! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
65!
66!
67! L.Fairhead 02/2000
68!
69! input:
70!   itime        numero du pas de temps
71!   klon         nombre total de points de grille
72!   dtime        pas de temps de la physique (en s)
73!   jour         jour dans l'annee en cours
74!   nisurf       index de la surface a traiter (1 = sol continental)
75!   knon         nombre de points de la surface a traiter
76!   knindex      index des points de la surface a traiter
77!   rlon         longitudes
78!   rlat         latitudes
79!   debut        logical: 1er appel a la physique
80!   lafin        logical: dernier appel a la physique
81!   ok_veget     logical: appel ou non au schema de surface continental
82!                     (si false calcul simplifie des fluxs sur les continents)
83!   zlev         hauteur de la premiere couche
84!   zlflu        epaisseur de la premier couche
85!   u1_lay       vitesse u 1ere couche
86!   v1_lay       vitesse v 1ere couche
87!   temp_air     temperature de l'air 1ere couche
88!   spechum      humidite specifique 1ere couche
89!   hum_air      humidite de l'air
90!   ccanopy      concentration CO2 canopee
91!   tq_cdrag     cdrag
92!   petAcoef     coeff. A de la resolution de la CL pour t
93!   peqAcoef     coeff. A de la resolution de la CL pour q
94!   petBcoef     coeff. B de la resolution de la CL pour t
95!   peqBcoef     coeff. B de la resolution de la CL pour q
96!   precip_rain  precipitation liquide
97!   precip_snow  precipitation solide
98!   lwdown       flux IR entrant a la surface
99!   swnet        flux solaire net
100!   swdown       flux solaire entrant a la surface
101!   tsurf        temperature de surface
102!   p1lay        pression 1er niveau (milieu de couche)
103!   cal          capacite calorifique du sol
104!   beta         evap reelle
105!   coef1lay     coefficient d'echange
106!   ps           pression au sol
107!   radsol       rayonnement net aus sol (LW + SW)
108!   dif_grnd     coeff. diffusion vers le sol profond
109!   ocean        type d'ocean utilise (force, slab, couple)
110!
111! output:
112!   evap         evaporation totale
113!   fluxsens     flux de chaleur sensible
114!   fluxlat      flux de chaleur latente
115!   tsol_rad     
116!   tsurf_new    temperature au sol
117!   alb_new      albedo
118!   emis_new     emissivite
119!   z0_new       surface roughness
120
121
122! Parametres d'entree
123  integer, intent(IN) :: itime
124  integer, intent(IN) :: klon
125  real, intent(IN) :: dtime
126  integer, intent(IN) :: jour
127  integer, intent(IN) :: nisurf
128  integer, intent(IN) :: knon
129  integer, dimension(knon), intent(in) :: knindex
130  logical, intent(IN) :: debut, lafin, ok_veget
131  real, dimension(klon), intent(IN) :: rlon, rlat
132  real, dimension(knon), intent(IN) :: zlev, zlflu
133  real, dimension(knon), intent(IN) :: u1_lay, v1_lay
134  real, dimension(knon), intent(IN) :: temp_air, spechum
135  real, dimension(knon), intent(IN) :: hum_air, ccanopy
136  real, dimension(knon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
137  real, dimension(knon), intent(IN) :: petBcoef, peqBcoef
138  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
139  real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps
140  real, dimension(knon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
141  real, dimension(knon), intent(IN) :: radsol, dif_grnd
142  character (len = 6)  :: ocean
143
144! Parametres de sortie
145  real, dimension(knon), intent(OUT):: evap, fluxsens, fluxlat
146  real, dimension(knon), intent(OUT):: tsol_rad, tsurf_new, alb_new
147  real, dimension(knon), intent(OUT):: emis_new, z0_new
148  real, dimension(knon), intent(OUT):: dflux_l, dflux_s
149  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
150
151! Local
152  character (len = 20) :: modname = 'interfsurf_hq'
153  character (len = 80) :: abort_message
154  logical, save        :: first_call = .true.
155  integer              :: error
156  logical              :: check = .true.
157
158  if (check) write(*,*) 'Entree ', modname
159!
160! On doit commencer par appeler les schemas de surfaces continentales
161! car l'ocean a besoin du ruissellement qui est y calcule
162!
163  if (first_call) then
164    if (nisurf /= is_ter .and. klon > 1) then
165      write(*,*)' *** Warning ***'
166      write(*,*)' nisurf = ',nisurf,' /= is_ter = ',is_ter
167      write(*,*)'or on doit commencer par les surfaces continentales'
168      abort_message='voir ci-dessus'
169      call abort_gcm(modname,abort_message,1)
170    endif
171    if (ocean /= 'slab  ' .and. ocean /= 'force ' .and. ocean /= 'couple') then
172      write(*,*)' *** Warning ***'
173      write(*,*)'Option couplage pour l''ocean = ', ocean
174      abort_message='option pour l''ocean non valable'
175      call abort_gcm(modname,abort_message,1)
176    endif
177  endif
178  first_call = .false.
179
180! Aiguillage vers les differents schemas de surface
181
182  if (nisurf == is_ter) then
183!
184! Surface "terre" appel a l'interface avec les sols continentaux
185!
186! allocation du run-off
187    if (.not. allocated(run_off)) then
188      allocate(run_off(knon), stat = error)
189      if (error /= 0) then
190        abort_message='Pb allocation run_off'
191        call abort_gcm(modname,abort_message,1)
192      endif
193    else if (size(run_off) /= knon) then
194      write(*,*)'Bizarre, le nombre de points continentaux'
195      write(*,*)'a change entre deux appels. Je continue ...'
196      deallocate(run_off, stat = error)
197      allocate(run_off(knon), stat = error)
198      if (error /= 0) then
199        abort_message='Pb allocation run_off'
200        call abort_gcm(modname,abort_message,1)
201      endif
202    endif
203!
204    call interfsol(itime, klon, dtime, nisurf, knon, &
205     &  knindex, rlon, rlat, &
206     &  debut, lafin, ok_veget, &
207     &  zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, &
208     &  tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
209     &  precip_rain, precip_snow, lwdown, swnet, swdown, &
210     &  tsurf, p1lay, cal, beta, coef1lay, ps, radsol, dif_grnd, &
211     &  evap, fluxsens, fluxlat, &             
212     &  tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
213     
214!
215  else if (nisurf == is_oce) then
216
217    if (check) write(*,*)'ocean, nisurf = ',nisurf
218!
219! Surface "ocean" appel a l'interface avec l'ocean
220!
221!    if (ocean == 'couple') then
222!      call interfoce(nisurf, ocean)
223!    else if (ocean == 'slab  ') then
224!      call interfoce(nisurf)
225!    else                              ! lecture conditions limites
226!      call interfoce(itime, dtime, jour, &
227!     &  klon, nisurf, knon, knindex, &
228!     &  debut, &
229!     &  tsurf_new, alb_new, z0_new, pctsrf_new)
230!    endif
231!
232    call calcul_fluxs( knon, dtime, &
233     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
234     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
235     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
236     &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
237!
238  else if (nisurf == is_sic) then
239
240    if (check) write(*,*)'sea ice, nisurf = ',nisurf
241
242!
243! Surface "glace de mer" appel a l'interface avec l'ocean
244!
245!    call interfoce(nisurf, ocean)
246!
247    call calcul_fluxs( knon, dtime, &
248     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
249     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
250     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
251     &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
252
253  else if (nisurf == is_lic) then
254
255    if (check) write(*,*)'glacier, nisurf = ',nisurf
256
257!
258! Surface "glacier continentaux" appel a l'interface avec le sol
259!
260!    call interfsol(nisurf)
261    call calcul_fluxs( knon, dtime, &
262     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
263     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
264     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
265     &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
266
267  else
268    write(*,*)'Index surface = ',nisurf
269    abort_message = 'Index surface non valable'
270    call abort_gcm(modname,abort_message,1)
271  endif
272
273  END SUBROUTINE interfsurf_hq
274
275!
276!#########################################################################
277!
278  SUBROUTINE interfsurf_vent(nisurf, knon   &         
279  &                     )
280!
281! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general
282! (sols continentaux, oceans, glaces) pour les tensions de vents.
283! En pratique l'interface se fait entre la couche limite du modele
284! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
285!
286!
287! L.Fairhead 02/2000
288!
289! input:
290!   nisurf       index de la surface a traiter (1 = sol continental)
291!   knon         nombre de points de la surface a traiter
292
293! Parametres d'entree
294  integer, intent(IN) :: nisurf
295  integer, intent(IN) :: knon
296
297
298  return
299  END SUBROUTINE interfsurf_vent
300!
301!#########################################################################
302!
303  SUBROUTINE interfsol(itime, klon, dtime, nisurf, knon, &
304     & knindex, rlon, rlat, &
305     & debut, lafin, ok_veget, &
306     & zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, &
307     & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
308     & precip_rain, precip_snow, lwdown, swnet, swdown, &
309     & tsurf, p1lay, cal, beta, coef1lay, ps, radsol, dif_grnd, &
310     & evap, fluxsens, fluxlat, &             
311     & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
312
313! Cette routine sert d'interface entre le modele atmospherique et le
314! modele de sol continental. Appel a sechiba
315!
316! L. Fairhead 02/2000
317!
318! input:
319!   itime        numero du pas de temps
320!   klon         nombre total de points de grille
321!   dtime        pas de temps de la physique (en s)
322!   nisurf       index de la surface a traiter (1 = sol continental)
323!   knon         nombre de points de la surface a traiter
324!   knindex      index des points de la surface a traiter
325!   rlon         longitudes de la grille entiere
326!   rlat         latitudes de la grille entiere
327!   debut        logical: 1er appel a la physique (lire les restart)
328!   lafin        logical: dernier appel a la physique (ecrire les restart)
329!   ok_veget     logical: appel ou non au schema de surface continental
330!                     (si false calcul simplifie des fluxs sur les continents)
331!   zlev         hauteur de la premiere couche
332!   zlflu       
333!   u1_lay       vitesse u 1ere couche
334!   v1_lay       vitesse v 1ere couche
335!   temp_air     temperature de l'air 1ere couche
336!   spechum      humidite specifique 1ere couche
337!   hum_air      humidite de l'air
338!   ccanopy      concentration CO2 canopee
339!   tq_cdrag     cdrag
340!   petAcoef     coeff. A de la resolution de la CL pour t
341!   peqAcoef     coeff. A de la resolution de la CL pour q
342!   petBcoef     coeff. B de la resolution de la CL pour t
343!   peqBcoef     coeff. B de la resolution de la CL pour q
344!   precip_rain  precipitation liquide
345!   precip_snow  precipitation solide
346!   lwdown       flux IR entrant a la surface
347!   swnet        flux solaire net
348!   swdown       flux solaire entrant a la surface
349!   tsurf        temperature de surface
350!   p1lay        pression 1er niveau (milieu de couche)
351!   cal          capacite calorifique du sol
352!   beta         evap reelle
353!   coef1lay     coefficient d'echange
354!   ps           pression au sol
355!   radsol       rayonnement net aus sol (LW + SW)
356!   dif_grnd     coeff. diffusion vers le sol profond     
357!   
358!
359! input/output
360!   run_off      ruissellement total
361!
362! output:
363!   evap         evaporation totale
364!   fluxsens     flux de chaleur sensible
365!   fluxlat      flux de chaleur latente
366!   tsol_rad     
367!   tsurf_new    temperature au sol
368!   alb_new      albedo
369!   emis_new     emissivite
370!   z0_new       surface roughness
371
372! Parametres d'entree
373  integer, intent(IN) :: itime
374  integer, intent(IN) :: klon
375  real, intent(IN)    :: dtime
376  integer, intent(IN) :: nisurf
377  integer, intent(IN) :: knon
378  integer, dimension(knon), intent(IN) :: knindex
379  logical, intent(IN) :: debut, lafin, ok_veget
380  real, dimension(klon), intent(IN) :: rlon, rlat
381  real, dimension(knon), intent(IN) :: zlev, zlflu
382  real, dimension(knon), intent(IN) :: u1_lay, v1_lay
383  real, dimension(knon), intent(IN) :: temp_air, spechum
384  real, dimension(knon), intent(IN) :: hum_air, ccanopy
385  real, dimension(knon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
386  real, dimension(knon), intent(IN) :: petBcoef, peqBcoef
387  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
388  real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps
389  real, dimension(knon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
390  real, dimension(knon), intent(IN) :: radsol, dif_grnd
391! Parametres de sortie
392  real, dimension(knon), intent(OUT):: evap, fluxsens, fluxlat
393  real, dimension(knon), intent(OUT):: tsol_rad, tsurf_new, alb_new
394  real, dimension(knon), intent(OUT):: emis_new, z0_new
395  real, dimension(knon), intent(OUT):: dflux_s, dflux_l
396
397! Local
398!
399  integer              :: ii
400  integer              :: error
401  character (len = 20) :: modname = 'interfsol'
402  character (len = 80) :: abort_message
403  logical              :: check = .true.
404! type de couplage dans sechiba
405!  character (len=10)   :: coupling = 'implicit'
406! drapeaux controlant les appels dans SECHIBA
407!  type(control_type), save   :: control_in
408! coordonnees geographiques
409  real, allocatable, dimension(:,:), save :: lalo
410! pts voisins
411  integer,allocatable, dimension(:,:), save :: neighbours
412! resolution de la grille
413  real, allocatable, dimension (:,:), save :: resolution
414! Identifieurs des fichiers restart et histoire
415  integer, save          :: rest_id, hist_id
416  integer, save          :: rest_id_stom, hist_id_stom
417
418  if (check) write(*,*)'Entree ', modname
419  if (check) write(*,*)'ok_veget = ',ok_veget
420
421! initialisation
422  if (.not. ok_veget) then
423    call calcul_fluxs( knon, dtime, &
424     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
425     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
426     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
427     &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
428  else
429!    if (debut) then
430!      !
431!      ! Configuration de parametres specifiques a la SSL
432!      !
433!      call intsurf_config(control_in)
434!      !
435!      ! Allouer et initialiser le tableau de coordonnees du sol
436!      !
437!      if (( .not. allocated(lalo))) then
438!        allocate(lalo(knon,2), stat = error)
439!        if (error /= 0) then
440!          abort_message='Pb allocation lalo'
441!          call abort_gcm(modname,abort_message,1)
442!        endif     
443!      endif
444!      do ii = 1, knon
445!        lalo(ii,1) = rlat(knindex(ii))
446!        lalo(ii,2) = rlon(knindex(ii))
447!      enddo
448      !-
449      !- Compute variable to help describe the grid
450      !- once the points are gathered.
451      !-
452!      IF ( (.NOT.ALLOCATED(neighbours))) THEN
453!        ALLOCATE(neighbours(knon,4), stat = error)
454!        if (error /= 0) then
455!          abort_message='Pb allocation neighbours'
456!          call abort_gcm(modname,abort_message,1)
457!        endif
458!      ENDIF
459!      IF ( (.NOT.ALLOCATED(resolution))) THEN
460!        ALLOCATE(resolution(knon,2), stat = error)
461!        if (error /= 0) then
462!          abort_message='Pb allocation resolution'
463!          call abort_gcm(modname,abort_message,1)
464!        endif
465!      ENDIF
466
467! call grid_stuff
468! call sechiba_restart_init
469! call sechiba_history_init
470
471!    endif                          ! (fin debut)
472
473!
474! Appel a la routine sols continentaux
475!
476
477!    call sechiba_main(itime, klon, knon, knindex, dtime, &
478!     & debut, lafin, coupling, control_in, &
479!     & lalo, neighbours, resolution,&
480!     & zlev, zlflu, u1_lay, v1_lay, spechum, temp_air,hum_air , ccanopy, &
481!     & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
482!     & precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
483!     & evap, fluxsens, fluxlat, &
484!     & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, &
485!     & rest_id, hist_id, rest_id_stom, hist_id_stom)
486
487!
488! Sauvegarde dans fichiers histoire
489!
490  endif  ! fin ok_veget
491
492  END SUBROUTINE interfsol
493!
494!#########################################################################
495!
496  SUBROUTINE interfoce_cpl(nisurf, ocean)
497
498! Cette routine sert d'interface entre le modele atmospherique et un
499! coupleur avec un modele d'ocean 'complet'
500!
501! L. Fairhead 02/2000
502!
503! input:
504!   nisurf       index de la surface a traiter (1 = sol continental)
505!   ocean        type d'ocean
506!
507! output:
508!
509
510! Parametres d'entree
511  integer, intent(IN) :: nisurf
512  character (len = 6)  :: ocean
513
514! Parametres de sortie
515
516! Variables locales
517
518
519! Initialisation
520! fichier restart et fichiers histoires
521
522! calcul des fluxs a passer
523
524  END SUBROUTINE interfoce_cpl
525!
526!#########################################################################
527!
528  SUBROUTINE interfoce_slab(nisurf)
529
530! Cette routine sert d'interface entre le modele atmospherique et un
531! modele de 'slab' ocean
532!
533! L. Fairhead 02/2000
534!
535! input:
536!   nisurf       index de la surface a traiter (1 = sol continental)
537!
538! output:
539!
540
541! Parametres d'entree
542  integer, intent(IN) :: nisurf
543
544  END SUBROUTINE interfoce_slab
545!
546!#########################################################################
547!
548  SUBROUTINE interfoce_lim(itime, dtime, jour, &
549     & klon, nisurf, knon, knindex, &
550     & debut,  &
551     & lmt_sst, lmt_alb, lmt_rug, pctsrf_new)
552
553! Cette routine sert d'interface entre le modele atmospherique et un fichier
554! de conditions aux limites
555!
556! L. Fairhead 02/2000
557!
558! input:
559!   itime        numero du pas de temps courant
560!   dtime        pas de temps de la physique (en s)
561!   jour         jour a lire dans l'annee
562!   nisurf       index de la surface a traiter (1 = sol continental)
563!   knon         nombre de points dans le domaine a traiter
564!   knindex      index des points de la surface a traiter
565!   klon         taille de la grille
566!   debut        logical: 1er appel a la physique (initialisation)
567!
568! output:
569!   lmt_sst      SST lues dans le fichier de CL
570!   lmt_alb      Albedo lu
571!   lmt_rug      longueur de rugosité lue
572!   pctsrf_new   sous-maille fractionnelle
573!
574
575
576! Parametres d'entree
577  integer, intent(IN) :: itime
578  real   , intent(IN) :: dtime
579  integer, intent(IN) :: jour
580  integer, intent(IN) :: nisurf
581  integer, intent(IN) :: knon
582  integer, intent(IN) :: klon
583  integer, dimension(knon), intent(in) :: knindex
584  logical, intent(IN) :: debut
585
586! Parametres de sortie
587  real, intent(out), dimension(knon) :: lmt_sst
588  real, intent(out), dimension(knon) :: lmt_alb
589  real, intent(out), dimension(knon) :: lmt_rug
590  real, intent(out), dimension(klon,nbsrf) :: pctsrf_new
591
592! Variables locales
593  integer     :: ii
594  integer     :: lmt_pas     ! frequence de lecture des conditions limites
595                             ! (en pas de physique)
596  logical,save :: deja_lu    ! pour indiquer que le jour a lire a deja
597                             ! lu pour une surface precedente
598  integer,save :: jour_lu
599  integer      :: ierr
600  character (len = 20) :: modname = 'interfoce_lim'
601  character (len = 80) :: abort_message
602  character (len = 20) :: fich ='limit'
603  logical     :: newlmt = .false.
604  logical     :: check = .true.
605! Champs lus dans le fichier de CL
606  real, allocatable , save, dimension(:) :: sst_lu, alb_lu, rug_lu, nat_lu
607  real, allocatable , save, dimension(:,:) :: pct_tmp
608!
609! quelques variables pour netcdf
610!
611#include "netcdf.inc"
612  integer              :: nid, nvarid
613  integer, dimension(2) :: start, epais
614!
615! Fin déclaration
616!
617   
618  if (debut) then
619    lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
620    jour_lu = jour - 1
621    allocate(sst_lu(klon))
622    allocate(alb_lu(klon))
623    allocate(rug_lu(klon))
624    allocate(nat_lu(klon))
625    allocate(pct_tmp(klon,nbsrf))
626  endif
627
628  if ((jour - jour_lu) /= 0) deja_lu = .false.
629 
630  if (check) write(*,*)modname,' :: jour_lu, deja_lu', jour_lu, deja_lu
631
632! Tester d'abord si c'est le moment de lire le fichier
633  if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then
634!
635! Ouverture du fichier
636!
637    ierr = NF_OPEN (fich, NF_NOWRITE,nid)
638    if (ierr.NE.NF_NOERR) then
639      abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
640      call abort_gcm(modname,abort_message,1)
641    endif
642!
643! La tranche de donnees a lire:
644!
645    start(1) = 1
646    start(2) = jour + 1
647    epais(1) = klon
648    epais(2) = 1
649!
650    if (newlmt) then
651!
652! Fraction "ocean"
653!
654      ierr = NF_INQ_VARID(nid, 'FOCE', nvarid)
655      if (ierr /= NF_NOERR) then
656        abort_message = 'Le champ <FOCE> est absent'
657        call abort_gcm(modname,abort_message,1)
658      endif
659#ifdef NC_DOUBLE
660      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pctsrf_new(1,is_oce))
661#else
662      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pctsrf_new(1,is_oce))
663#endif
664      if (ierr /= NF_NOERR) then
665        abort_message = 'Lecture echouee pour <FOCE>'
666        call abort_gcm(modname,abort_message,1)
667      endif
668!
669! Fraction "glace de mer"
670!
671      ierr = NF_INQ_VARID(nid, 'FSIC', nvarid)
672      if (ierr /= NF_NOERR) then
673        abort_message = 'Le champ <FSIC> est absent'
674        call abort_gcm(modname,abort_message,1)
675      endif
676#ifdef NC_DOUBLE
677      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pctsrf_new(1,is_sic))
678#else
679      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pctsrf_new(1,is_sic))
680#endif
681      if (ierr /= NF_NOERR) then
682        abort_message = 'Lecture echouee pour <FSIC>'
683        call abort_gcm(modname,abort_message,1)
684      endif
685!
686! Fraction "terre"
687!
688      ierr = NF_INQ_VARID(nid, 'FTER', nvarid)
689      if (ierr /= NF_NOERR) then
690        abort_message = 'Le champ <FTER> est absent'
691        call abort_gcm(modname,abort_message,1)
692      endif
693#ifdef NC_DOUBLE
694      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pctsrf_new(1,is_ter))
695#else
696      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pctsrf_new(1,is_ter))
697#endif
698      if (ierr /= NF_NOERR) then
699        abort_message = 'Lecture echouee pour <FTER>'
700        call abort_gcm(modname,abort_message,1)
701      endif
702!
703! Fraction "glacier terre"
704!
705      ierr = NF_INQ_VARID(nid, 'FLIC', nvarid)
706      if (ierr /= NF_NOERR) then
707        abort_message = 'Le champ <FLIC> est absent'
708        call abort_gcm(modname,abort_message,1)
709      endif
710#ifdef NC_DOUBLE
711      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pctsrf_new(1,is_lic))
712#else
713      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pctsrf_new(1,is_lic))
714#endif
715      if (ierr /= NF_NOERR) then
716        abort_message = 'Lecture echouee pour <FLIC>'
717        call abort_gcm(modname,abort_message,1)
718      endif
719!
720    else  ! on en est toujours a rnatur
721!
722      ierr = NF_INQ_VARID(nid, 'NAT', nvarid)
723      if (ierr /= NF_NOERR) then
724        abort_message = 'Le champ <NAT> est absent'
725        call abort_gcm(modname,abort_message,1)
726      endif
727#ifdef NC_DOUBLE
728      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, nat_lu)
729#else
730      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, nat_lu)
731#endif
732      if (ierr /= NF_NOERR) then
733        abort_message = 'Lecture echouee pour <NAT>'
734        call abort_gcm(modname,abort_message,1)
735      endif
736!
737! Remplissage des fractions de surface
738! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice
739!
740      pct_tmp = 0.0
741      do ii = 1, klon
742        pct_tmp(ii,nint(nat_lu(ii)) + 1) = 1.
743      enddo
744
745!
746!  On se retrouve avec ocean en 1 et terre en 2 alors qu'on veut le contraire
747!
748      pctsrf_new = pct_tmp
749      pctsrf_new (:,2)= pct_tmp (:,1)
750      pctsrf_new (:,1)= pct_tmp (:,2)
751      pct_tmp = pctsrf_new
752    endif ! fin test sur newlmt
753!
754! Lecture SST
755!
756    ierr = NF_INQ_VARID(nid, 'SST', nvarid)
757    if (ierr /= NF_NOERR) then
758      abort_message = 'Le champ <SST> est absent'
759      call abort_gcm(modname,abort_message,1)
760    endif
761#ifdef NC_DOUBLE
762    ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, sst_lu)
763#else
764    ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, sst_lu)
765#endif
766    if (ierr /= NF_NOERR) then
767      abort_message = 'Lecture echouee pour <SST>'
768      call abort_gcm(modname,abort_message,1)
769    endif   
770!
771! Lecture Albedo
772!
773    ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
774    if (ierr /= NF_NOERR) then
775      abort_message = 'Le champ <ALB> est absent'
776      call abort_gcm(modname,abort_message,1)
777    endif
778#ifdef NC_DOUBLE
779    ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, alb_lu)
780#else
781    ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, alb_lu)
782#endif
783    if (ierr /= NF_NOERR) then
784      abort_message = 'Lecture echouee pour <ALB>'
785      call abort_gcm(modname,abort_message,1)
786    endif
787!
788! Lecture rugosité
789!
790    ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
791    if (ierr /= NF_NOERR) then
792      abort_message = 'Le champ <RUG> est absent'
793      call abort_gcm(modname,abort_message,1)
794    endif
795#ifdef NC_DOUBLE
796    ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, rug_lu)
797#else
798    ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)
799#endif
800    if (ierr /= NF_NOERR) then
801      abort_message = 'Lecture echouee pour <RUG>'
802      call abort_gcm(modname,abort_message,1)
803    endif
804
805!
806! Fin de lecture
807!
808    ierr = NF_CLOSE(nid)
809    deja_lu = .true.
810    jour_lu = jour
811  endif
812!
813! Recopie des variables dans les champs de sortie
814!
815  do ii = 1, knon
816    lmt_sst(ii) = sst_lu(knindex(ii))
817    lmt_alb(ii) = alb_lu(knindex(ii))
818    lmt_rug(ii) = rug_lu(knindex(ii))
819  enddo
820  pctsrf_new = pct_tmp
821
822  END SUBROUTINE interfoce_lim
823
824!
825!#########################################################################
826!
827
828  SUBROUTINE calcul_fluxs( knon, dtime, &
829     & tsurf, p1lay, cal, beta, coef1lay, ps, &
830     & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
831     & petAcoef, peqAcoef, petBcoef, peqBcoef, &
832     & tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
833
834! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
835! une temperature de surface (au cas ou ok_veget = false)
836!
837! L. Fairhead 4/2000
838!
839! input:
840!   knon         nombre de points a traiter
841!   tsurf        temperature de surface
842!   p1lay        pression 1er niveau (milieu de couche)
843!   cal          capacite calorifique du sol
844!   beta         evap reelle
845!   coef1lay     coefficient d'echange
846!   ps           pression au sol
847!   petAcoef     coeff. A de la resolution de la CL pour t
848!   peqAcoef     coeff. A de la resolution de la CL pour q
849!   petBcoef     coeff. B de la resolution de la CL pour t
850!   peqBcoef     coeff. B de la resolution de la CL pour q
851!   radsol       rayonnement net aus sol (LW + SW)
852!   dif_grnd     coeff. diffusion vers le sol profond
853!
854! output:
855!   tsurf_new    temperature au sol
856!   fluxsens     flux de chaleur sensible
857!   fluxlat      flux de chaleur latente
858!   dflux_s      derivee du flux de chaleur sensible / Ts
859!   dflux_l      derivee du flux de chaleur latente  / Ts
860!
861
862#include "YOMCST.inc"
863#include "YOETHF.inc"
864#include "FCTTRE.inc"
865
866! Parametres d'entree
867  integer, intent(IN) :: knon
868  real   , intent(IN) :: dtime
869  real, dimension(knon), intent(IN) :: petAcoef, peqAcoef
870  real, dimension(knon), intent(IN) :: petBcoef, peqBcoef
871  real, dimension(knon), intent(IN) :: ps, q1lay
872  real, dimension(knon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
873  real, dimension(knon), intent(IN) :: radsol, dif_grnd
874  real, dimension(knon), intent(IN) :: t1lay, u1lay, v1lay
875
876! Parametres sorties
877  real, dimension(knon), intent(OUT):: tsurf_new, fluxlat, fluxsens
878  real, dimension(knon), intent(OUT):: dflux_s, dflux_l
879
880! Variables locales
881  integer :: i
882  real, dimension(knon) :: zx_mh, zx_nh, zx_oh
883  real, dimension(knon) :: zx_mq, zx_nq, zx_oq
884  real, dimension(knon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
885  real, dimension(knon) :: zx_sl, zx_k1,  zx_dq,  zx_cq,  zx_dh, zx_ch
886  real, dimension(knon) :: zx_h_ts, zx_q_0 , d_ts
887  real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
888  real, parameter :: t_grnd = 271.35, t_coup = 273.15
889  logical         :: check = .true.
890  character (len = 20)  :: modname = 'calcul_fluxs'
891
892  if (check) write(*,*)'Entree ', modname
893!
894! Initialisation
895!
896!
897! zx_qs = qsat en kg/kg
898!
899  DO i = 1, knon
900    zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
901  IF (thermcep) THEN
902      zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
903      zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
904      zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
905      zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
906      zx_qs=MIN(0.5,zx_qs)
907      zcor=1./(1.-retv*zx_qs)
908      zx_qs=zx_qs*zcor
909      zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
910     &                 /RLVTT / zx_pkh(i)
911    ELSE
912      IF (tsurf(i).LT.t_coup) THEN
913        zx_qs = qsats(tsurf(i)) / ps(i)
914        zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
915     &                    / zx_pkh(i)
916      ELSE
917        zx_qs = qsatl(tsurf(i)) / ps(i)
918        zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
919     &               / zx_pkh(i)
920      ENDIF
921    ENDIF
922    zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
923    zx_qsat(i) = zx_qs
924    zx_coef(i) = coef1lay(i) &
925     & * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
926     & * p1lay(i)/(RD*t1lay(i))
927   
928  ENDDO
929
930
931! === Calcul de la temperature de surface ===
932!
933! zx_sl = chaleur latente d'evaporation ou de sublimation
934!
935  do i = 1, knon
936    zx_sl(i) = RLVTT
937    if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
938    zx_k1(i) = zx_coef(i)
939  enddo
940
941
942  do i = 1, knon
943! Q
944    zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
945    zx_mq(i) = beta(i) * zx_k1(i) * &
946     &             (peqAcoef(i) - zx_qsat(i) &
947     &                          + zx_dq_s_dt(i) * tsurf(i)) &
948     &             / zx_oq(i)
949    zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
950     &                              / zx_oq(i)
951
952! H
953    zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
954    zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
955    zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
956
957! Tsurface
958    tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
959     &             (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
960     &                 + dif_grnd(i) * t_grnd * dtime)/ &
961     &          ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
962     &                       zx_nh(i) + zx_sl(i) * zx_nq(i)) & 
963     &                     + dtime * dif_grnd(i))
964    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
965    d_ts(i) = tsurf_new(i) - tsurf(i)
966    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
967!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
968!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
969    fluxlat(i) = zx_mq(i) + zx_nq(i) * tsurf_new(i)
970    fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
971! Derives des flux dF/dTs (W m-2 K-1):
972    dflux_s(i) = zx_nh(i)
973    dflux_l(i) = (zx_sl(i) * zx_nq(i))
974  ENDDO
975
976  END SUBROUTINE calcul_fluxs
977!
978!#########################################################################
979!
980
981  END MODULE interface_surf
Note: See TracBrowser for help on using the repository browser.