source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyaqua_mod.F90 @ 3960

Last change on this file since 3960 was 3960, checked in by ymipsl, 8 years ago

Fix uninitialized variables in LMDZ physics.

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