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

Last change on this file since 4080 was 4033, checked in by millour, 8 years ago

Add some missing initializations for the aquaplanet case.
EM
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
339!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
340!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
341!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
342    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
343    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
344
345!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
346!ym probably the uninitialized value was 0 for standard (regular grid) case
347    falb_dif(:,:,:)=0
348
349
350    CALL phyredem('startphy.nc')
351
352    PRINT *, 'iniaqua: after phyredem'
353    CALL phys_state_var_end
354
355    RETURN
356  END SUBROUTINE iniaqua
357
358
359  ! ====================================================================
360  ! ====================================================================
361  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
362    USE dimphy
363    IMPLICIT NONE
364    ! ====================================================================
365    ! =============================================================
366    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
367    ! Auteur : A. Campoy et F. Hourdin
368    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
369    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
370    ! connaissant la declinaison, la latitude et la longitude.
371
372    ! Dans cette version particuliere, on calcule le rayonnement
373    ! moyen sur l'année à chaque latitude.
374    ! angle zenithal calculé pour obtenir un
375    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
376    ! en moyenne annuelle.
377    ! Spécifique de la terre. Utilisé pour les aqua planetes.
378
379    ! Rque   : Different de la routine angle en ce sens que zenang
380    ! fournit des moyennes de pmu0 et non des valeurs
381    ! instantanees, du coup frac prend toutes les valeurs
382    ! entre 0 et 1.
383    ! Date   : premiere version le 13 decembre 1994
384    ! revu pour  GCM  le 30 septembre 1996
385    ! ===============================================================
386    ! longi----INPUT : la longitude vraie de la terre dans son plan
387    ! solaire a partir de l'equinoxe de printemps (degre)
388    ! gmtime---INPUT : temps universel en fraction de jour
389    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
390    ! lat------INPUT : latitude en degres
391    ! long-----INPUT : longitude en degres
392    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
393    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
394    ! ================================================================
395    include "YOMCST.h"
396    ! ================================================================
397    LOGICAL cycle_diurne
398    REAL gmtime
399    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
400    ! ================================================================
401    INTEGER i
402    REAL gmtime1, gmtime2
403    REAL pi_local
404
405
406    REAL rmu0m(klon), rmu0a(klon)
407
408
409    pi_local = 4.0*atan(1.0)
410
411    ! ================================================================
412    ! Calcul de l'angle zenithal moyen sur la journee
413    ! ================================================================
414
415    DO i = 1, klon
416      fract(i) = 1.
417      ! Calcule du flux moyen
418      IF (abs(rlat(i))<=28.75) THEN
419        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
420      ELSE IF (abs(rlat(i))<=43.75) THEN
421        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
422      ELSE IF (abs(rlat(i))<=71.25) THEN
423        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
424      ELSE
425        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
426      END IF
427    END DO
428
429    ! ================================================================
430    ! Avec ou sans cycle diurne
431    ! ================================================================
432
433    IF (cycle_diurne) THEN
434
435      ! On redecompose flux  au sommet suivant un cycle diurne idealise
436      ! identique a toutes les latitudes.
437
438      DO i = 1, klon
439        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
440        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
441          rmu0a(i)/sqrt(2.)
442      END DO
443
444      DO i = 1, klon
445        IF (rmu0(i)<=0.) THEN
446          rmu0(i) = 0.
447          fract(i) = 0.
448        ELSE
449          fract(i) = 1.
450        END IF
451      END DO
452
453      ! Affichage de l'angel zenitale
454      ! print*,'************************************'
455      ! print*,'************************************'
456      ! print*,'************************************'
457      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
458      ! print*,'rmu0m=',rmu0m(i)
459      ! print*,'rmu0a=',rmu0a(i)
460      ! print*,'rmu0=',rmu0(i)
461
462    ELSE
463
464      DO i = 1, klon
465        fract(i) = 0.5
466        rmu0(i) = rmu0m(i)*2.
467      END DO
468
469    END IF
470
471    RETURN
472  END SUBROUTINE zenang_an
473
474  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
475
476  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
477      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
478
479    USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
480    USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
481#ifdef CPP_XIOS
482    USE xios
483#endif
484    IMPLICIT NONE
485
486    include "netcdf.inc"
487
488    INTEGER, INTENT (IN) :: klon
489    REAL, INTENT (IN) :: phy_nat(klon, 360)
490    REAL, INTENT (IN) :: phy_alb(klon, 360)
491    REAL, INTENT (IN) :: phy_sst(klon, 360)
492    REAL, INTENT (IN) :: phy_bil(klon, 360)
493    REAL, INTENT (IN) :: phy_rug(klon, 360)
494    REAL, INTENT (IN) :: phy_ice(klon, 360)
495    REAL, INTENT (IN) :: phy_fter(klon, 360)
496    REAL, INTENT (IN) :: phy_foce(klon, 360)
497    REAL, INTENT (IN) :: phy_flic(klon, 360)
498    REAL, INTENT (IN) :: phy_fsic(klon, 360)
499
500    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
501      ! on the whole physics grid
502 
503#ifdef CPP_XIOS
504    PRINT *, 'writelim: Ecriture du fichier limit'
505
506    CALL gather_omp(phy_foce, phy_mpi)
507    IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
508
509    CALL gather_omp(phy_fsic, phy_mpi)
510    IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
511     
512    CALL gather_omp(phy_fter, phy_mpi)
513    IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
514     
515    CALL gather_omp(phy_flic, phy_mpi)
516    IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
517
518    CALL gather_omp(phy_sst, phy_mpi)
519    IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
520
521    CALL gather_omp(phy_bil, phy_mpi)
522    IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
523
524    CALL gather_omp(phy_alb, phy_mpi)
525    IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
526
527    CALL gather_omp(phy_rug, phy_mpi)
528    IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
529#endif
530  END SUBROUTINE writelim_unstruct
531
532
533
534  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
535      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
536
537    USE mod_phys_lmdz_para, ONLY: is_master
538    USE mod_grid_phy_lmdz, ONLY: klon_glo
539    USE mod_phys_lmdz_transfert_para, ONLY: gather
540    IMPLICIT NONE
541    include "netcdf.inc"
542
543    INTEGER, INTENT (IN) :: klon
544    REAL, INTENT (IN) :: phy_nat(klon, 360)
545    REAL, INTENT (IN) :: phy_alb(klon, 360)
546    REAL, INTENT (IN) :: phy_sst(klon, 360)
547    REAL, INTENT (IN) :: phy_bil(klon, 360)
548    REAL, INTENT (IN) :: phy_rug(klon, 360)
549    REAL, INTENT (IN) :: phy_ice(klon, 360)
550    REAL, INTENT (IN) :: phy_fter(klon, 360)
551    REAL, INTENT (IN) :: phy_foce(klon, 360)
552    REAL, INTENT (IN) :: phy_flic(klon, 360)
553    REAL, INTENT (IN) :: phy_fsic(klon, 360)
554
555    REAL :: phy_glo(klon_glo, 360) ! temporary variable, to store phy_***(:)
556      ! on the whole physics grid
557    INTEGER :: k
558    INTEGER ierr
559    INTEGER dimfirst(3)
560    INTEGER dimlast(3)
561
562    INTEGER nid, ndim, ntim
563    INTEGER dims(2), debut(2), epais(2)
564    INTEGER id_tim
565    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
566    INTEGER id_fter, id_foce, id_fsic, id_flic
567
568    IF (is_master) THEN
569
570      PRINT *, 'writelim: Ecriture du fichier limit'
571
572      ierr = nf_create('limit.nc', NF_NETCDF4, nid)
573
574      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
575        'Fichier conditions aux limites')
576      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
577      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
578      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
579
580      dims(1) = ndim
581      dims(2) = ntim
582
583#ifdef NC_DOUBLE
584      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
585#else
586      ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
587#endif
588      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
589
590#ifdef NC_DOUBLE
591      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
592#else
593      ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
594#endif
595      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
596        'Nature du sol (0,1,2,3)')
597
598#ifdef NC_DOUBLE
599      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
600#else
601      ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
602#endif
603      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
604        'Temperature superficielle de la mer')
605
606#ifdef NC_DOUBLE
607      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
608#else
609      ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
610#endif
611      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
612        'Reference flux de chaleur au sol')
613
614#ifdef NC_DOUBLE
615      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
616#else
617      ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
618#endif
619      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
620
621#ifdef NC_DOUBLE
622      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
623#else
624      ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
625#endif
626      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
627
628#ifdef NC_DOUBLE
629      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
630#else
631      ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
632#endif
633      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
634#ifdef NC_DOUBLE
635      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
636#else
637      ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
638#endif
639      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
640#ifdef NC_DOUBLE
641      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
642#else
643      ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
644#endif
645      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
646#ifdef NC_DOUBLE
647      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
648#else
649      ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
650#endif
651      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
652
653      ierr = nf_enddef(nid)
654      IF (ierr/=nf_noerr) THEN
655        WRITE (*, *) 'writelim error: failed to end define mode'
656        WRITE (*, *) nf_strerror(ierr)
657      END IF
658
659
660      ! write the 'times'
661      DO k = 1, 360
662#ifdef NC_DOUBLE
663        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
664#else
665        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
666#endif
667        IF (ierr/=nf_noerr) THEN
668          WRITE (*, *) 'writelim error with temps(k),k=', k
669          WRITE (*, *) nf_strerror(ierr)
670        END IF
671      END DO
672
673    END IF ! of if (is_master)
674
675    ! write the fields, after having collected them on master
676
677    CALL gather(phy_nat, phy_glo)
678    IF (is_master) THEN
679#ifdef NC_DOUBLE
680      ierr = nf_put_var_double(nid, id_nat, phy_glo)
681#else
682      ierr = nf_put_var_real(nid, id_nat, phy_glo)
683#endif
684      IF (ierr/=nf_noerr) THEN
685        WRITE (*, *) 'writelim error with phy_nat'
686        WRITE (*, *) nf_strerror(ierr)
687      END IF
688    END IF
689
690    CALL gather(phy_sst, phy_glo)
691    IF (is_master) THEN
692#ifdef NC_DOUBLE
693      ierr = nf_put_var_double(nid, id_sst, phy_glo)
694#else
695      ierr = nf_put_var_real(nid, id_sst, phy_glo)
696#endif
697      IF (ierr/=nf_noerr) THEN
698        WRITE (*, *) 'writelim error with phy_sst'
699        WRITE (*, *) nf_strerror(ierr)
700      END IF
701    END IF
702
703    CALL gather(phy_bil, phy_glo)
704    IF (is_master) THEN
705#ifdef NC_DOUBLE
706      ierr = nf_put_var_double(nid, id_bils, phy_glo)
707#else
708      ierr = nf_put_var_real(nid, id_bils, phy_glo)
709#endif
710      IF (ierr/=nf_noerr) THEN
711        WRITE (*, *) 'writelim error with phy_bil'
712        WRITE (*, *) nf_strerror(ierr)
713      END IF
714    END IF
715
716    CALL gather(phy_alb, phy_glo)
717    IF (is_master) THEN
718#ifdef NC_DOUBLE
719      ierr = nf_put_var_double(nid, id_alb, phy_glo)
720#else
721      ierr = nf_put_var_real(nid, id_alb, phy_glo)
722#endif
723      IF (ierr/=nf_noerr) THEN
724        WRITE (*, *) 'writelim error with phy_alb'
725        WRITE (*, *) nf_strerror(ierr)
726      END IF
727    END IF
728
729    CALL gather(phy_rug, phy_glo)
730    IF (is_master) THEN
731#ifdef NC_DOUBLE
732      ierr = nf_put_var_double(nid, id_rug, phy_glo)
733#else
734      ierr = nf_put_var_real(nid, id_rug, phy_glo)
735#endif
736      IF (ierr/=nf_noerr) THEN
737        WRITE (*, *) 'writelim error with phy_rug'
738        WRITE (*, *) nf_strerror(ierr)
739      END IF
740    END IF
741
742    CALL gather(phy_fter, phy_glo)
743    IF (is_master) THEN
744#ifdef NC_DOUBLE
745      ierr = nf_put_var_double(nid, id_fter, phy_glo)
746#else
747      ierr = nf_put_var_real(nid, id_fter, phy_glo)
748#endif
749      IF (ierr/=nf_noerr) THEN
750        WRITE (*, *) 'writelim error with phy_fter'
751        WRITE (*, *) nf_strerror(ierr)
752      END IF
753    END IF
754
755    CALL gather(phy_foce, phy_glo)
756    IF (is_master) THEN
757#ifdef NC_DOUBLE
758      ierr = nf_put_var_double(nid, id_foce, phy_glo)
759#else
760      ierr = nf_put_var_real(nid, id_foce, phy_glo)
761#endif
762      IF (ierr/=nf_noerr) THEN
763        WRITE (*, *) 'writelim error with phy_foce'
764        WRITE (*, *) nf_strerror(ierr)
765      END IF
766    END IF
767
768    CALL gather(phy_fsic, phy_glo)
769    IF (is_master) THEN
770#ifdef NC_DOUBLE
771      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
772#else
773      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
774#endif
775      IF (ierr/=nf_noerr) THEN
776        WRITE (*, *) 'writelim error with phy_fsic'
777        WRITE (*, *) nf_strerror(ierr)
778      END IF
779    END IF
780
781    CALL gather(phy_flic, phy_glo)
782    IF (is_master) THEN
783#ifdef NC_DOUBLE
784      ierr = nf_put_var_double(nid, id_flic, phy_glo)
785#else
786      ierr = nf_put_var_real(nid, id_flic, phy_glo)
787#endif
788      IF (ierr/=nf_noerr) THEN
789        WRITE (*, *) 'writelim error with phy_flic'
790        WRITE (*, *) nf_strerror(ierr)
791      END IF
792    END IF
793
794    ! close file:
795    IF (is_master) THEN
796      ierr = nf_close(nid)
797    END IF
798
799  END SUBROUTINE writelim
800
801  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
802
803  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
804    USE dimphy
805    IMPLICIT NONE
806
807    INTEGER nlon, type_profil, i, k, j
808    REAL :: rlatd(nlon), phy_sst(nlon, 360)
809    INTEGER imn, imx, amn, amx, kmn, kmx
810    INTEGER p, pplus, nlat_max
811    PARAMETER (nlat_max=72)
812    REAL x_anom_sst(nlat_max)
813
814    IF (klon/=nlon) STOP 'probleme de dimensions dans iniaqua'
815    WRITE (*, *) ' profil_sst: type_profil=', type_profil
816    DO i = 1, 360
817      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
818
819      ! Rajout fbrlmd
820
821      IF (type_profil==1) THEN
822        ! Méthode 1 "Control" faible plateau à l'Equateur
823        DO j = 1, klon
824          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
825          ! PI/3=1.047197551
826          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
827            phy_sst(j, i) = 273.
828          END IF
829        END DO
830      END IF
831      IF (type_profil==2) THEN
832        ! Méthode 2 "Flat" fort plateau à l'Equateur
833        DO j = 1, klon
834          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
835          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
836            phy_sst(j, i) = 273.
837          END IF
838        END DO
839      END IF
840
841
842      IF (type_profil==3) THEN
843        ! Méthode 3 "Qobs" plateau réel à l'Equateur
844        DO j = 1, klon
845          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
846            rlatd(j))**4)
847          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
848            phy_sst(j, i) = 273.
849          END IF
850        END DO
851      END IF
852
853      IF (type_profil==4) THEN
854        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
855        DO j = 1, klon
856          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
857            rlatd(j))**4)
858          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
859            phy_sst(j, i) = 273.
860          END IF
861        END DO
862      END IF
863
864      IF (type_profil==5) THEN
865        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
866        DO j = 1, klon
867          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
868            *rlatd(j))**4)
869          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
870            phy_sst(j, i) = 273. + 2.
871          END IF
872
873        END DO
874      END IF
875
876      IF (type_profil==6) THEN
877        ! Méthode 6 "cst" valeur constante de SST
878        DO j = 1, klon
879          phy_sst(j, i) = 288.
880        END DO
881      END IF
882
883
884      IF (type_profil==7) THEN
885        ! Méthode 7 "cst" valeur constante de SST +2
886        DO j = 1, klon
887          phy_sst(j, i) = 288. + 2.
888        END DO
889      END IF
890
891      p = 0
892      IF (type_profil==8) THEN
893        ! Méthode 8 profil anomalies SST du modèle couplé AR4
894        DO j = 1, klon
895          IF (rlatd(j)==rlatd(j-1)) THEN
896            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
897              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
898            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
899              phy_sst(j, i) = 273. + x_anom_sst(pplus)
900            END IF
901          ELSE
902            p = p + 1
903            pplus = 73 - p
904            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
905              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
906            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
907              phy_sst(j, i) = 273. + x_anom_sst(pplus)
908            END IF
909            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
910          END IF
911        END DO
912      END IF
913
914      IF (type_profil==9) THEN
915        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
916        DO j = 1, klon
917          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
918            *rlatd(j))**4)
919          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
920            phy_sst(j, i) = 273. - 2.
921          END IF
922        END DO
923      END IF
924
925
926      IF (type_profil==10) THEN
927        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
928        DO j = 1, klon
929          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
930            *rlatd(j))**4)
931          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
932            phy_sst(j, i) = 273. + 4.
933          END IF
934        END DO
935      END IF
936
937      IF (type_profil==11) THEN
938        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
939        DO j = 1, klon
940          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
941            rlatd(j))**4)
942          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
943            phy_sst(j, i) = 273.
944          END IF
945        END DO
946      END IF
947
948      IF (type_profil==12) THEN
949        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
950        DO j = 1, klon
951          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
952            *rlatd(j))**4)
953          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
954            phy_sst(j, i) = 273. + 4.
955          END IF
956        END DO
957      END IF
958
959      IF (type_profil==13) THEN
960        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
961        DO j = 1, klon
962          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
963            rlatd(j))**4)
964          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
965            phy_sst(j, i) = 273.
966          END IF
967        END DO
968      END IF
969
970      IF (type_profil==14) THEN
971        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
972        DO j = 1, klon
973          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
974            *rlatd(j))**4)
975          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
976            phy_sst(j, i) = 273.
977          END IF
978        END DO
979      END IF
980
981      if (type_profil.EQ.20) then
982      print*,'Profile SST 20'
983!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
984
985      do j=1,klon
986        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
987      enddo
988      endif
989
990      if (type_profil.EQ.21) then
991      print*,'Profile SST 21'
992!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
993      do j=1,klon
994        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
995      enddo
996      endif
997
998
999
1000    END DO
1001
1002    ! IM beg : verif profil SST: phy_sst
1003    amn = min(phy_sst(1,1), 1000.)
1004    amx = max(phy_sst(1,1), -1000.)
1005    imn = 1
1006    kmn = 1
1007    imx = 1
1008    kmx = 1
1009    DO k = 1, 360
1010      DO i = 2, nlon
1011        IF (phy_sst(i,k)<amn) THEN
1012          amn = phy_sst(i, k)
1013          imn = i
1014          kmn = k
1015        END IF
1016        IF (phy_sst(i,k)>amx) THEN
1017          amx = phy_sst(i, k)
1018          imx = i
1019          kmx = k
1020        END IF
1021      END DO
1022    END DO
1023
1024    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
1025    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
1026    ! IM end : verif profil SST: phy_sst
1027
1028    RETURN
1029  END SUBROUTINE profil_sst
1030
1031END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.