source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phylmd/phyaqua_mod.F90 @ 3603

Last change on this file since 3603 was 3312, checked in by Laurent Fairhead, 7 years ago

Continuing phasing of DYNAMICO and LMDZ physics
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 30.5 KB
Line 
1MODULE phyaqua_mod
2  ! Routines complementaires pour la physique planetaire.
3  IMPLICIT NONE
4
5CONTAINS
6
7  SUBROUTINE iniaqua(nlon, iflag_phys)
8
9    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10    ! Creation d'un etat initial et de conditions aux limites
11    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
12    ! du modele LMDZ dans sa version terrestre.
13    ! iflag_phys est un parametre qui controle
14    ! iflag_phys = N
15    ! de 100 a 199 : aqua planetes avec SST forcees
16    ! N-100 determine le type de SSTs
17    ! de 200 a 299 : terra planetes avec Ts calcule
18
19    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20
21    USE dimphy, ONLY: klon
22    USE geometry_mod, ONLY : latitude
23    USE surface_data, ONLY: type_ocean, ok_veget
24    USE pbl_surface_mod, ONLY: pbl_surface_init
25    USE fonte_neige_mod, ONLY: fonte_neige_init
26    USE phys_state_var_mod
27    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
28                                day_ini,day_end
29    USE indice_sol_mod
30    USE nrtype, ONLY: pi
31!    USE ioipsl
32    USE mod_phys_lmdz_para, ONLY: is_master
33    USE mod_phys_lmdz_transfert_para, ONLY: bcast
34    USE mod_grid_phy_lmdz
35    USE ioipsl_getin_p_mod, ONLY : getin_p
36    IMPLICIT NONE
37
38    include "YOMCST.h"
39    include "clesphys.h"
40    include "dimsoil.h"
41
42    INTEGER, INTENT (IN) :: nlon, iflag_phys
43    ! IM ajout latfi, lonfi
44!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
45
46    INTEGER type_profil, type_aqua
47
48    ! Ajouts initialisation des surfaces
49    REAL :: run_off_lic_0(nlon)
50    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
51    REAL :: tsoil(nlon, nsoilmx, nbsrf)
52    REAL :: tslab(nlon), seaice(nlon)
53    REAL fder(nlon)
54
55
56
57    ! Arguments :
58    ! -----------
59
60    ! integer radpas
61    INTEGER it, unit, i, k, itap
62
63    REAL rugos, albedo
64    REAL tsurf
65    REAL time, timestep, day, day0
66    REAL qsol_f
67    REAL rugsrel(nlon)
68    LOGICAL alb_ocean
69
70    CHARACTER *80 ans, file_forctl, file_fordat, file_start
71    CHARACTER *100 file, var
72    CHARACTER *2 cnbl
73
74    REAL phy_nat(nlon, 360)
75    REAL phy_alb(nlon, 360)
76    REAL phy_sst(nlon, 360)
77    REAL phy_bil(nlon, 360)
78    REAL phy_rug(nlon, 360)
79    REAL phy_ice(nlon, 360)
80    REAL phy_fter(nlon, 360)
81    REAL phy_foce(nlon, 360)
82    REAL phy_fsic(nlon, 360)
83    REAL phy_flic(nlon, 360)
84
85    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
86!$OMP THREADPRIVATE(read_climoz)
87
88    ! -------------------------------------------------------------------------
89    ! declaration pour l'appel a phyredem
90    ! -------------------------------------------------------------------------
91
92    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
93    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
94    ! real pbl_tke(nlon,llm,nbsrf)
95    ! real rain_fall(nlon),snow_fall(nlon)
96    ! real solsw(nlon), sollw(nlon),radsol(nlon)
97    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
98    ! real ratqs(nlon,llm)
99    ! real clwcon(nlon,llm)
100
101    INTEGER longcles
102    PARAMETER (longcles=20)
103    REAL clesphy0(longcles)
104
105
106    ! -----------------------------------------------------------------------
107    ! dynamial tendencies :
108    ! ---------------------
109
110    INTEGER l, ierr, aslun
111
112    REAL paire
113
114
115    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
116    ! INITIALISATIONS
117    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
118
119    ! -----------------------------------------------------------------------
120    ! Initialisations  des constantes
121    ! -------------------------------
122
123
124    type_aqua = iflag_phys/100
125    type_profil = iflag_phys - type_aqua*100
126    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
127
128    IF (klon/=nlon) THEN
129      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
130      STOP 'probleme de dimensions dans iniaqua'
131    END IF
132    CALL phys_state_var_init(read_climoz)
133
134
135    read_climoz = 0
136    day0 = 217.
137    day = day0
138    it = 0
139    time = 0.
140
141    ! -----------------------------------------------------------------------
142    ! initialisations de la physique
143    ! -----------------------------------------------------------------------
144
145    day_ini = day_ref
146    day_end = day_ini + ndays
147
148    nbapp_rad = 24
149    CALL getin_p('nbapp_rad', nbapp_rad)
150
151    ! ---------------------------------------------------------------------
152    ! Creation des conditions aux limites:
153    ! ------------------------------------
154    ! Initialisations des constantes
155    ! Ajouter les manquants dans planete.def... (albedo etc)
156    co2_ppm = 348.
157    CALL getin_p('co2_ppm', co2_ppm)
158
159    solaire = 1365.
160    CALL getin_p('solaire', solaire)
161 
162    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
163    ! type_aqua
164    alb_ocean = .TRUE.
165    CALL getin_p('alb_ocean', alb_ocean)
166
167    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
168    WRITE (*, *) 'iniaqua: solaire=', solaire
169    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
170
171    radsol = 0.
172    qsol_f = 10.
173
174    ! Conditions aux limites:
175    ! -----------------------
176
177    qsol(:) = qsol_f
178    rugsrel = 0.0 ! (rugsrel = rugoro)
179    rugoro = 0.0
180    u_ancien = 0.0
181    v_ancien = 0.0
182    agesno = 50.0
183    ! Relief plat
184    zmea = 0.
185    zstd = 0.
186    zsig = 0.
187    zgam = 0.
188    zthe = 0.
189    zpic = 0.
190    zval = 0.
191
192    ! Une seule surface
193    pctsrf = 0.
194    IF (type_aqua==1) THEN
195      rugos = 1.E-4
196      albedo = 0.19
197      pctsrf(:, is_oce) = 1.
198    ELSE IF (type_aqua==2) THEN
199      rugos = 0.03
200      albedo = 0.1
201      pctsrf(:, is_ter) = 1.
202    END IF
203
204    CALL getin_p('rugos', rugos)
205
206    WRITE (*, *) 'iniaqua: rugos=', rugos
207    zmasq(:) = pctsrf(:, is_ter)
208
209    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
210    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
211
212    ! Si alb_ocean on calcule un albedo oceanique moyen
213    ! if (alb_ocean) then
214    ! Voir pourquoi on avait ca.
215    ! CALL ini_alb_oce(phy_alb)
216    ! else
217    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
218    ! endif !alb_ocean
219
220    DO i = 1, 360
221      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
222      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
223
224      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
225      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
226      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
227      phy_ice(:, i) = 0.0 ! fraction de glace (?)
228      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
229      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
230      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
231      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
232    END DO
233    ! IM calcul profil sst
234    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
235
236    IF (grid_type==unstructured) THEN
237      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
238                             phy_fter, phy_foce, phy_flic, phy_fsic)
239    ELSE
240     
241       CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
242                     phy_fter, phy_foce, phy_flic, phy_fsic)
243    ENDIF
244
245    ! ---------------------------------------------------------------------
246    ! Ecriture de l'etat initial:
247    ! ---------------------------
248
249
250    ! Ecriture etat initial physique
251
252    timestep = pdtphys
253    radpas = nint(rday/timestep/float(nbapp_rad))
254
255    DO i = 1, longcles
256      clesphy0(i) = 0.
257    END DO
258    clesphy0(1) = float(iflag_con)
259    clesphy0(2) = float(nbapp_rad)
260    ! IF( cycle_diurne  ) clesphy0(3) =  1.
261    clesphy0(3) = 1. ! cycle_diurne
262    clesphy0(4) = 1. ! soil_model
263    clesphy0(5) = 1. ! new_oliq
264    clesphy0(6) = 0. ! ok_orodr
265    clesphy0(7) = 0. ! ok_orolf
266    clesphy0(8) = 0. ! ok_limitvrai
267
268
269    ! =======================================================================
270    ! Profils initiaux
271    ! =======================================================================
272
273    ! On initialise les temperatures de surfaces comme les sst
274    DO i = 1, nlon
275      ftsol(i, :) = phy_sst(i, 1)
276      tsoil(i, :, :) = phy_sst(i, 1)
277      tslab(i) = phy_sst(i, 1)
278    END DO
279
280    falbe(:, :) = albedo
281    falblw(:, :) = albedo
282    rain_fall(:) = 0.
283    snow_fall(:) = 0.
284    solsw(:) = 0.
285    sollw(:) = 0.
286    radsol(:) = 0.
287
288    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
289    ! intialisation bidon mais pas grave
290    t_ancien(:, :) = 0.
291    q_ancien(:, :) = 0.
292    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
293    rnebcon = 0.
294    ratqs = 0.
295    clwcon = 0.
296    pbl_tke = 1.E-8
297
298    ! variables supplementaires pour appel a plb_surface_init
299    fder(:) = 0.
300    seaice(:) = 0.
301    run_off_lic_0 = 0.
302    fevap = 0.
303
304
305    ! Initialisations necessaires avant phyredem
306    type_ocean = 'force'
307    CALL fonte_neige_init(run_off_lic_0)
308    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
309    snsrf(:, :) = 0. ! couverture de neige des sous surface
310    z0m(:, :) = rugos ! couverture de neige des sous surface
311    z0h=z0m
312
313
314    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
315
316    PRINT *, 'iniaqua: before phyredem'
317
318    pbl_tke(:,:,:) = 1.e-8
319    falb1 = albedo
320    falb2 = albedo
321    zmax0 = 0.
322    f0 = 0.
323    sig1 = 0.
324    w01 = 0.
325    wake_deltat = 0.
326    wake_deltaq = 0.
327    wake_s = 0.
328    wake_dens = 0.
329    wake_cstar = 0.
330    wake_pe = 0.
331    wake_fip = 0.
332    fm_therm = 0.
333    entr_therm = 0.
334    detr_therm = 0.
335    ale_bl = 0.
336    ale_bl_trig =0.
337    alp_bl =0.
338    treedrg(:,:,:)=0.
339
340!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
341!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
342!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
343    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
344    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
345
346!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
347!ym probably the uninitialized value was 0 for standard (regular grid) case
348    falb_dif(:,:,:)=0
349
350
351    CALL phyredem('startphy.nc')
352
353    PRINT *, 'iniaqua: after phyredem'
354    CALL phys_state_var_end
355
356    RETURN
357  END SUBROUTINE iniaqua
358
359
360  ! ====================================================================
361  ! ====================================================================
362  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
363    USE dimphy
364    IMPLICIT NONE
365    ! ====================================================================
366    ! =============================================================
367    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
368    ! Auteur : A. Campoy et F. Hourdin
369    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
370    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
371    ! connaissant la declinaison, la latitude et la longitude.
372
373    ! Dans cette version particuliere, on calcule le rayonnement
374    ! moyen sur l'année à chaque latitude.
375    ! angle zenithal calculé pour obtenir un
376    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
377    ! en moyenne annuelle.
378    ! Spécifique de la terre. Utilisé pour les aqua planetes.
379
380    ! Rque   : Different de la routine angle en ce sens que zenang
381    ! fournit des moyennes de pmu0 et non des valeurs
382    ! instantanees, du coup frac prend toutes les valeurs
383    ! entre 0 et 1.
384    ! Date   : premiere version le 13 decembre 1994
385    ! revu pour  GCM  le 30 septembre 1996
386    ! ===============================================================
387    ! longi----INPUT : la longitude vraie de la terre dans son plan
388    ! solaire a partir de l'equinoxe de printemps (degre)
389    ! gmtime---INPUT : temps universel en fraction de jour
390    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
391    ! lat------INPUT : latitude en degres
392    ! long-----INPUT : longitude en degres
393    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
394    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
395    ! ================================================================
396    include "YOMCST.h"
397    ! ================================================================
398    LOGICAL cycle_diurne
399    REAL gmtime
400    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
401    ! ================================================================
402    INTEGER i
403    REAL gmtime1, gmtime2
404    REAL pi_local
405
406
407    REAL rmu0m(klon), rmu0a(klon)
408
409
410    pi_local = 4.0*atan(1.0)
411
412    ! ================================================================
413    ! Calcul de l'angle zenithal moyen sur la journee
414    ! ================================================================
415
416    DO i = 1, klon
417      fract(i) = 1.
418      ! Calcule du flux moyen
419      IF (abs(rlat(i))<=28.75) THEN
420        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
421      ELSE IF (abs(rlat(i))<=43.75) THEN
422        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
423      ELSE IF (abs(rlat(i))<=71.25) THEN
424        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
425      ELSE
426        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
427      END IF
428    END DO
429
430    ! ================================================================
431    ! Avec ou sans cycle diurne
432    ! ================================================================
433
434    IF (cycle_diurne) THEN
435
436      ! On redecompose flux  au sommet suivant un cycle diurne idealise
437      ! identique a toutes les latitudes.
438
439      DO i = 1, klon
440        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
441        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
442          rmu0a(i)/sqrt(2.)
443      END DO
444
445      DO i = 1, klon
446        IF (rmu0(i)<=0.) THEN
447          rmu0(i) = 0.
448          fract(i) = 0.
449        ELSE
450          fract(i) = 1.
451        END IF
452      END DO
453
454      ! Affichage de l'angel zenitale
455      ! print*,'************************************'
456      ! print*,'************************************'
457      ! print*,'************************************'
458      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
459      ! print*,'rmu0m=',rmu0m(i)
460      ! print*,'rmu0a=',rmu0a(i)
461      ! print*,'rmu0=',rmu0(i)
462
463    ELSE
464
465      DO i = 1, klon
466        fract(i) = 0.5
467        rmu0(i) = rmu0m(i)*2.
468      END DO
469
470    END IF
471
472    RETURN
473  END SUBROUTINE zenang_an
474
475  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
476
477  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
478      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
479
480    USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
481    USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
482#ifdef CPP_XIOS
483    USE xios
484#endif
485    IMPLICIT NONE
486
487    include "netcdf.inc"
488
489    INTEGER, INTENT (IN) :: klon
490    REAL, INTENT (IN) :: phy_nat(klon, 360)
491    REAL, INTENT (IN) :: phy_alb(klon, 360)
492    REAL, INTENT (IN) :: phy_sst(klon, 360)
493    REAL, INTENT (IN) :: phy_bil(klon, 360)
494    REAL, INTENT (IN) :: phy_rug(klon, 360)
495    REAL, INTENT (IN) :: phy_ice(klon, 360)
496    REAL, INTENT (IN) :: phy_fter(klon, 360)
497    REAL, INTENT (IN) :: phy_foce(klon, 360)
498    REAL, INTENT (IN) :: phy_flic(klon, 360)
499    REAL, INTENT (IN) :: phy_fsic(klon, 360)
500
501    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
502      ! on the whole physics grid
503 
504#ifdef CPP_XIOS
505    PRINT *, 'writelim: Ecriture du fichier limit'
506
507    CALL gather_omp(phy_foce, phy_mpi)
508    IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
509
510    CALL gather_omp(phy_fsic, phy_mpi)
511    IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
512     
513    CALL gather_omp(phy_fter, phy_mpi)
514    IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
515     
516    CALL gather_omp(phy_flic, phy_mpi)
517    IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
518
519    CALL gather_omp(phy_sst, phy_mpi)
520    IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
521
522    CALL gather_omp(phy_bil, phy_mpi)
523    IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
524
525    CALL gather_omp(phy_alb, phy_mpi)
526    IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
527
528    CALL gather_omp(phy_rug, phy_mpi)
529    IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
530#endif
531  END SUBROUTINE writelim_unstruct
532
533
534
535  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
536      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
537
538    USE mod_phys_lmdz_para, ONLY: is_master
539    USE mod_grid_phy_lmdz, ONLY: klon_glo
540    USE mod_phys_lmdz_transfert_para, ONLY: gather
541    IMPLICIT NONE
542    include "netcdf.inc"
543
544    INTEGER, INTENT (IN) :: klon
545    REAL, INTENT (IN) :: phy_nat(klon, 360)
546    REAL, INTENT (IN) :: phy_alb(klon, 360)
547    REAL, INTENT (IN) :: phy_sst(klon, 360)
548    REAL, INTENT (IN) :: phy_bil(klon, 360)
549    REAL, INTENT (IN) :: phy_rug(klon, 360)
550    REAL, INTENT (IN) :: phy_ice(klon, 360)
551    REAL, INTENT (IN) :: phy_fter(klon, 360)
552    REAL, INTENT (IN) :: phy_foce(klon, 360)
553    REAL, INTENT (IN) :: phy_flic(klon, 360)
554    REAL, INTENT (IN) :: phy_fsic(klon, 360)
555
556    REAL :: phy_glo(klon_glo, 360) ! temporary variable, to store phy_***(:)
557      ! on the whole physics grid
558    INTEGER :: k
559    INTEGER ierr
560    INTEGER dimfirst(3)
561    INTEGER dimlast(3)
562
563    INTEGER nid, ndim, ntim
564    INTEGER dims(2), debut(2), epais(2)
565    INTEGER id_tim
566    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
567    INTEGER id_fter, id_foce, id_fsic, id_flic
568
569    IF (is_master) THEN
570
571      PRINT *, 'writelim: Ecriture du fichier limit'
572
573      ierr = nf_create('limit.nc', NF_NETCDF4, nid)
574
575      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
576        'Fichier conditions aux limites')
577      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
578      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
579      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
580
581      dims(1) = ndim
582      dims(2) = ntim
583
584#ifdef NC_DOUBLE
585      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
586#else
587      ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
588#endif
589      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
590
591#ifdef NC_DOUBLE
592      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
593#else
594      ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
595#endif
596      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
597        'Nature du sol (0,1,2,3)')
598
599#ifdef NC_DOUBLE
600      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
601#else
602      ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
603#endif
604      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
605        'Temperature superficielle de la mer')
606
607#ifdef NC_DOUBLE
608      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
609#else
610      ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
611#endif
612      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
613        'Reference flux de chaleur au sol')
614
615#ifdef NC_DOUBLE
616      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
617#else
618      ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
619#endif
620      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
621
622#ifdef NC_DOUBLE
623      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
624#else
625      ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
626#endif
627      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
628
629#ifdef NC_DOUBLE
630      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
631#else
632      ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
633#endif
634      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
635#ifdef NC_DOUBLE
636      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
637#else
638      ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
639#endif
640      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
641#ifdef NC_DOUBLE
642      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
643#else
644      ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
645#endif
646      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
647#ifdef NC_DOUBLE
648      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
649#else
650      ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
651#endif
652      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
653
654      ierr = nf_enddef(nid)
655      IF (ierr/=nf_noerr) THEN
656        WRITE (*, *) 'writelim error: failed to end define mode'
657        WRITE (*, *) nf_strerror(ierr)
658      END IF
659
660
661      ! write the 'times'
662      DO k = 1, 360
663#ifdef NC_DOUBLE
664        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
665#else
666        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
667#endif
668        IF (ierr/=nf_noerr) THEN
669          WRITE (*, *) 'writelim error with temps(k),k=', k
670          WRITE (*, *) nf_strerror(ierr)
671        END IF
672      END DO
673
674    END IF ! of if (is_master)
675
676    ! write the fields, after having collected them on master
677
678    CALL gather(phy_nat, phy_glo)
679    IF (is_master) THEN
680#ifdef NC_DOUBLE
681      ierr = nf_put_var_double(nid, id_nat, phy_glo)
682#else
683      ierr = nf_put_var_real(nid, id_nat, phy_glo)
684#endif
685      IF (ierr/=nf_noerr) THEN
686        WRITE (*, *) 'writelim error with phy_nat'
687        WRITE (*, *) nf_strerror(ierr)
688      END IF
689    END IF
690
691    CALL gather(phy_sst, phy_glo)
692    IF (is_master) THEN
693#ifdef NC_DOUBLE
694      ierr = nf_put_var_double(nid, id_sst, phy_glo)
695#else
696      ierr = nf_put_var_real(nid, id_sst, phy_glo)
697#endif
698      IF (ierr/=nf_noerr) THEN
699        WRITE (*, *) 'writelim error with phy_sst'
700        WRITE (*, *) nf_strerror(ierr)
701      END IF
702    END IF
703
704    CALL gather(phy_bil, phy_glo)
705    IF (is_master) THEN
706#ifdef NC_DOUBLE
707      ierr = nf_put_var_double(nid, id_bils, phy_glo)
708#else
709      ierr = nf_put_var_real(nid, id_bils, phy_glo)
710#endif
711      IF (ierr/=nf_noerr) THEN
712        WRITE (*, *) 'writelim error with phy_bil'
713        WRITE (*, *) nf_strerror(ierr)
714      END IF
715    END IF
716
717    CALL gather(phy_alb, phy_glo)
718    IF (is_master) THEN
719#ifdef NC_DOUBLE
720      ierr = nf_put_var_double(nid, id_alb, phy_glo)
721#else
722      ierr = nf_put_var_real(nid, id_alb, phy_glo)
723#endif
724      IF (ierr/=nf_noerr) THEN
725        WRITE (*, *) 'writelim error with phy_alb'
726        WRITE (*, *) nf_strerror(ierr)
727      END IF
728    END IF
729
730    CALL gather(phy_rug, phy_glo)
731    IF (is_master) THEN
732#ifdef NC_DOUBLE
733      ierr = nf_put_var_double(nid, id_rug, phy_glo)
734#else
735      ierr = nf_put_var_real(nid, id_rug, phy_glo)
736#endif
737      IF (ierr/=nf_noerr) THEN
738        WRITE (*, *) 'writelim error with phy_rug'
739        WRITE (*, *) nf_strerror(ierr)
740      END IF
741    END IF
742
743    CALL gather(phy_fter, phy_glo)
744    IF (is_master) THEN
745#ifdef NC_DOUBLE
746      ierr = nf_put_var_double(nid, id_fter, phy_glo)
747#else
748      ierr = nf_put_var_real(nid, id_fter, phy_glo)
749#endif
750      IF (ierr/=nf_noerr) THEN
751        WRITE (*, *) 'writelim error with phy_fter'
752        WRITE (*, *) nf_strerror(ierr)
753      END IF
754    END IF
755
756    CALL gather(phy_foce, phy_glo)
757    IF (is_master) THEN
758#ifdef NC_DOUBLE
759      ierr = nf_put_var_double(nid, id_foce, phy_glo)
760#else
761      ierr = nf_put_var_real(nid, id_foce, phy_glo)
762#endif
763      IF (ierr/=nf_noerr) THEN
764        WRITE (*, *) 'writelim error with phy_foce'
765        WRITE (*, *) nf_strerror(ierr)
766      END IF
767    END IF
768
769    CALL gather(phy_fsic, phy_glo)
770    IF (is_master) THEN
771#ifdef NC_DOUBLE
772      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
773#else
774      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
775#endif
776      IF (ierr/=nf_noerr) THEN
777        WRITE (*, *) 'writelim error with phy_fsic'
778        WRITE (*, *) nf_strerror(ierr)
779      END IF
780    END IF
781
782    CALL gather(phy_flic, phy_glo)
783    IF (is_master) THEN
784#ifdef NC_DOUBLE
785      ierr = nf_put_var_double(nid, id_flic, phy_glo)
786#else
787      ierr = nf_put_var_real(nid, id_flic, phy_glo)
788#endif
789      IF (ierr/=nf_noerr) THEN
790        WRITE (*, *) 'writelim error with phy_flic'
791        WRITE (*, *) nf_strerror(ierr)
792      END IF
793    END IF
794
795    ! close file:
796    IF (is_master) THEN
797      ierr = nf_close(nid)
798    END IF
799
800  END SUBROUTINE writelim
801
802  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
803
804  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
805    USE dimphy
806    IMPLICIT NONE
807
808    INTEGER nlon, type_profil, i, k, j
809    REAL :: rlatd(nlon), phy_sst(nlon, 360)
810    INTEGER imn, imx, amn, amx, kmn, kmx
811    INTEGER p, pplus, nlat_max
812    PARAMETER (nlat_max=72)
813    REAL x_anom_sst(nlat_max)
814
815    IF (klon/=nlon) STOP 'probleme de dimensions dans iniaqua'
816    WRITE (*, *) ' profil_sst: type_profil=', type_profil
817    DO i = 1, 360
818      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
819
820      ! Rajout fbrlmd
821
822      IF (type_profil==1) THEN
823        ! Méthode 1 "Control" faible plateau à l'Equateur
824        DO j = 1, klon
825          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
826          ! PI/3=1.047197551
827          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
828            phy_sst(j, i) = 273.
829          END IF
830        END DO
831      END IF
832      IF (type_profil==2) THEN
833        ! Méthode 2 "Flat" fort plateau à l'Equateur
834        DO j = 1, klon
835          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
836          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
837            phy_sst(j, i) = 273.
838          END IF
839        END DO
840      END IF
841
842
843      IF (type_profil==3) THEN
844        ! Méthode 3 "Qobs" plateau réel à l'Equateur
845        DO j = 1, klon
846          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
847            rlatd(j))**4)
848          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
849            phy_sst(j, i) = 273.
850          END IF
851        END DO
852      END IF
853
854      IF (type_profil==4) THEN
855        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
856        DO j = 1, klon
857          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
858            rlatd(j))**4)
859          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
860            phy_sst(j, i) = 273.
861          END IF
862        END DO
863      END IF
864
865      IF (type_profil==5) THEN
866        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
867        DO j = 1, klon
868          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
869            *rlatd(j))**4)
870          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
871            phy_sst(j, i) = 273. + 2.
872          END IF
873
874        END DO
875      END IF
876
877      IF (type_profil==6) THEN
878        ! Méthode 6 "cst" valeur constante de SST
879        DO j = 1, klon
880          phy_sst(j, i) = 288.
881        END DO
882      END IF
883
884
885      IF (type_profil==7) THEN
886        ! Méthode 7 "cst" valeur constante de SST +2
887        DO j = 1, klon
888          phy_sst(j, i) = 288. + 2.
889        END DO
890      END IF
891
892      p = 0
893      IF (type_profil==8) THEN
894        ! Méthode 8 profil anomalies SST du modèle couplé AR4
895        DO j = 1, klon
896          IF (rlatd(j)==rlatd(j-1)) THEN
897            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
898              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
899            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
900              phy_sst(j, i) = 273. + x_anom_sst(pplus)
901            END IF
902          ELSE
903            p = p + 1
904            pplus = 73 - p
905            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
906              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
907            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
908              phy_sst(j, i) = 273. + x_anom_sst(pplus)
909            END IF
910            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
911          END IF
912        END DO
913      END IF
914
915      IF (type_profil==9) THEN
916        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
917        DO j = 1, klon
918          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
919            *rlatd(j))**4)
920          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
921            phy_sst(j, i) = 273. - 2.
922          END IF
923        END DO
924      END IF
925
926
927      IF (type_profil==10) THEN
928        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
929        DO j = 1, klon
930          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
931            *rlatd(j))**4)
932          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
933            phy_sst(j, i) = 273. + 4.
934          END IF
935        END DO
936      END IF
937
938      IF (type_profil==11) THEN
939        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
940        DO j = 1, klon
941          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
942            rlatd(j))**4)
943          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
944            phy_sst(j, i) = 273.
945          END IF
946        END DO
947      END IF
948
949      IF (type_profil==12) THEN
950        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
951        DO j = 1, klon
952          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
953            *rlatd(j))**4)
954          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
955            phy_sst(j, i) = 273. + 4.
956          END IF
957        END DO
958      END IF
959
960      IF (type_profil==13) THEN
961        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
962        DO j = 1, klon
963          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
964            rlatd(j))**4)
965          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
966            phy_sst(j, i) = 273.
967          END IF
968        END DO
969      END IF
970
971      IF (type_profil==14) THEN
972        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
973        DO j = 1, klon
974          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
975            *rlatd(j))**4)
976          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
977            phy_sst(j, i) = 273.
978          END IF
979        END DO
980      END IF
981
982      if (type_profil.EQ.20) then
983      print*,'Profile SST 20'
984!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
985
986      do j=1,klon
987        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
988      enddo
989      endif
990
991      if (type_profil.EQ.21) then
992      print*,'Profile SST 21'
993!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
994      do j=1,klon
995        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
996      enddo
997      endif
998
999
1000
1001    END DO
1002
1003    ! IM beg : verif profil SST: phy_sst
1004    amn = min(phy_sst(1,1), 1000.)
1005    amx = max(phy_sst(1,1), -1000.)
1006    imn = 1
1007    kmn = 1
1008    imx = 1
1009    kmx = 1
1010    DO k = 1, 360
1011      DO i = 2, nlon
1012        IF (phy_sst(i,k)<amn) THEN
1013          amn = phy_sst(i, k)
1014          imn = i
1015          kmn = k
1016        END IF
1017        IF (phy_sst(i,k)>amx) THEN
1018          amx = phy_sst(i, k)
1019          imx = i
1020          kmx = k
1021        END IF
1022      END DO
1023    END DO
1024
1025    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
1026    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
1027    ! IM end : verif profil SST: phy_sst
1028
1029    RETURN
1030  END SUBROUTINE profil_sst
1031
1032END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.